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1.   INTRODUCTION:   VARIABLE  MULTISTEP  AND  MULTIVALUE  METHODS 

Consider  the  initial  value  problem  for  the  single  differential 
equation 

y'  =  f(y,t),  y(0)  =  yQ  (l.l) 

We  assume   that   f(y,t)    satisfies   appropriate  Lipschitz  and  continuity 
conditions   such  that   a  unique   solution  y(t)   exists   and  y(t)    e   C'[o,b]. 
Later  on   it  is   assumed  that  y(t)   has   as  many  higher  derivatives   as 
necessary. 

There   are  various  methods  to   compute   the  numerical   approximation 
of  y(t).     Among  them  are  the  general  k-step  methods   in  the  form 

k 

I      (a.    y      .    +  hg.    f      .)   =   0  (1.2) 

.    -        l     n-i  i     n-i 

i=0 

where  h  (a  constant)  E  step  size 

y   .  E  the  numerical  approximation  of  y(t)  at  t   .  e  (n-i)h 
n-i  n-i 

f  ..  E  f(y   .  ,  t   .) 
n-i      n-i   n-i 

a.    and  B.    are   constants, 
l  l 

It  has   long  been  established  that  by  specifying  an   appropriate 

set  of  values   for  a.    and  8.    one   can  obtain  a  stable   and  convergent  k-step 

11 

method.   Moreover,  there  exists  more  than  one  k-step  method.   For  instance, 
Adams-Bashforth  methods  (see  Henrici  [9],  Section  5.1),  Adams -Moul ton 
methods  (see  Henrici  [9],  Section  5.1),  and  stiffly  stable  methods 
(see  Gear  [8],  Section  11. l)  are  all  special  cases  of  (1.2). 

In  this  thesis  we  are  concerned  with  variable  step  size  mutlistep 
methods  which  are  generalizations  of  (1.2).   It  turns  out  that  a.  and  B. 


are  no  longer  constants.   Even  though  they  are  still  independent  of 


y   .  and  f   . ,  they  are  dependent  on  the  step  sizes  taken. 

We  are  also  going  to  consider  multivalue  methods  (see  Gear  [8], 
Section  T»l)  which  in  the  case  of  constant  step  size  h  have  the 
following  form 


4,(0)  ::  BV-i 


(1.3) 


Y  .     .»  •  Y  M  +  C  G(Y  /  %)  m  -  0.  1,...(1.M 
-n,(m+l)   -n,(m)   —   -n,(m)       •  *    ••*■'■•  ' 


where  B  is  a  constant  matrix,  C  a  constant  column  vector,  Y  ,  is  the 

—  -n-1 

column  vector  of  saved  information 

iT 


ln_!  »  tyn.r  yn- 


2»" 


yn.k»  hyn-r 


-  hyn-k] 


and 


G(Y     (    N)   =  h(f(y      f    ,,   t    )   -  y'    ,    J 
-n,(m)  wn,(m)'      n  n,(m) 

(The  description  of  this  procedure    (l.3)-(l.5)  will  be  elaborated  in 

Section  1.1. ) 

Predictor-corrector  multistep  methods   are   a  particular   case   of 

multivalue  methods.      For  instance,    a  U-value  Adams -Bashforth-Moul ton 

method  has  the   form 


(1.5) 


B  = 


1  23/12 

0  3 

0  1 

0  0 


-16/12  5/12 

-3  1 

0  0 

1  0 


C  r 


3/8 

1 
0 
0 


When  step  changing  is  necessary,  we   shall  consider  an  equivalent 
representation  of   (1.3)    and   (l.U)   in  an   attempt  to  make   step   changing 
simple.      This   is  based  on  the  use   of  the  Nordsieck  vector    [ll] 


a 
-n 


.1=  [Vr'Ci hk^l/k!lT 


instead  of  Y   .  Accordingly,  we  have  the  representation 


a   ,  v  =  A  a 
-n,(0)     -n-1 


(1.6) 

a   /  , . v  =  a   /  x  +  I    F(a   ,  J  (1.7) 

-n,(m+l)   -n,(m)   —   -n,(m) 

equivalent  to  (2.3)  and  (2.U).   Here,  again,  A  is  a  constant  matrix, 

(2)        (k) 

I   a  constant  vector,  F  the  same  function  as  G,  and  y'  n,  y,v  ' v 

-  "  '     Jn-1'  Jn-1  »•*•»  ^n-i 

are  numerical  derivatives  of  y  _  solely  used  in  the  multi value  method. 

n-1 

Note  that  for  instance  y1    %   f(y  _  ,  t  _  ) .   It  was  shown  by  Gear  [8] 
that  a  general  multistep  predictor-corrector  method  can  he  expressed 
by  representation  (1.6) -(1.7). 

To  demonstrate  the  correspondence  between  (l.3)-(lA)  and 

(l.6)-(l.7)  let  us  consider  U-value  Adams -Bashforth-Moulton  method.   In 

T 
this  case,  Y  .  is  [y  _,  hy!  .  ,  hy'  ...  hy'  _]  .   The  information  stored 
-n-1      n-1    n-1    n-2    n-3 

in  Y  ,  and  a  ,  determines  the  same  third  degree  interpolating  poly- 
-n-1     -n-1 


nomial  [  8  ]•   Therefore,  Y  .  and  a  .  are  related  by  a  transformati 

—n-1     -n-1 

T,  i.e. 


on 


a    =  T  Y  . 

-n-1     -n-1 


where 


T  = 


1 

0 

0 

0 

0 

1 

0 

0 

0 

3A 

-1 

iA 

0 

i/6 

-1/3 

1/6 

Consequently,  we  have  for  the  representation  (1.6) -(1.7) 


A  =   TBT 


,-1 


1111 
0  12  3 
0  0  13 
0     0     0      1 


3/8 

1 
I    =   TC  = 

3A 

1/6 


Thus,  our  approach  to  step  changing  is  two -fold.   On  the  one 
hand  we  look  at  the  interpolation  method  "based  on  multivalue  methods  in 
Nordsieck  form,  and  on  the  other  we  consider  a  variable  step  method  based 
on  (1.2).  We  shall  elaborate  on  these  two  approaches  here. 

1.1  Interpolation  Method  Using  Nordsieck  Vector 

Since  we  are  working  with  variable  step  size,  we  shall  assume 

throughout  this  paper  that  h    is  the  size  of  step  advancing  from  time 

t  .  to  t  ,  i.e.  h  -,=t  -t  _  .  We  adopt  the  conventional  use  of  y 
n-1     n        n-1    n    n-1  n-i 

to  indicate  the  numerical  approximation  of  the  solution  y(t)  at  t   .  and 

n-i 

naturally  f   .  =  f (y   . ,  t   . ) . 
n-i      n-i   n-i 

Step  changing  for  a  Nordsieck  vector  a    is  easily  accomplished 


by  the  multiplication  by  the  matrix 


'n-1 


■    o 

n-1 


n-1 


when  u 


n-1 


Formally , 


is  called  a  step  changing  factor, 

n-2 


(1.8) 


-n-1 


h 


yn-l 

n-1  yn-l 

• 
• 

= 

k     • 

n-1      (JO 

k!      ''n-l 



n-1 


o 


o 


n-1 


n-1 


h  n   y' 
n-2  °n-l 


h 


n-2   (k) 


k!   yn-l 


n-l-n-1 


(1.9) 


We  define  the  step  changing  algorithm  from  t   .  to  t   as  follows 

n-1  n 


1.      Prescribe   the   step   changing  factor   oo 


n-1 


2.   Use  (1.8)  to  obtain  C 


n-1' 


3.   Use  (1.9)  to  find  the  Nordsieck  vector  a*  . 

—n-1 


with  step  size  h 


n-1' 


h.      Replace  a     ,    in    (1.6)  by  a*   , . 
-n-1  —n-1 

Incorporating  this   algorithm  into    (l.6)-(l.7)  we  obtain 


^i,(0)   "  A  Cn-1  Vl 


a      ,_N    +  I    F(a_    /_^) 


-n,(m+l)        — n,(m)        -  1V2Ti,(m) 

This  modified  representation  (l.lO)-(l.ll )  is  called  an  interpolation 

method  since  the  step  changing  involves  an  interpolation  process. 

The  prediction  process  (1.10)  does  not  contain  any  knowledge 

of  the  differential  equation.   We  just  extrapolate  a   ,.x  from  a  _, 

— n,(,U;      -n-1 

we  need  to  correct  the  amount  by  which  a  , nS    does  not  satisfy  the 

-n,(0) 


(1.10) 

(1.11) 


So 


differential  equation;  this  is  performed  in  (l.ll).   Theoretically,  one 

may  repeat  the  correction  as  many  times  as  one  pleases  to  attain 

convergence.  However,  correction  steps  call  for  additional  computation 

and  about  two  such  steps  are  usually  recommended.   Suppose  a   ,  s    is  the 

value  at  the  final  correction  step,  then  a  x  a   ,,,>.. 

*'  -n   -n,(M) 

Suppose  the  matrix  A  is  a  Pascal  triangular  matrix  whose  (i,j) 
element  is  (  . )  ,  k  >  j  ■•_   i  >_  0,  and  the  step  size  is  kept  constant  through- 
out  the  numerical  integration,  then  there  exists  k-step  methods  of  maximum 
order  that  are  stable.   (For  constant  step  method,  the  meaning  of  order 
and  stability  are  in  accordance  with  Gear  [8]  or  Henrici  [9])-   A  proof 
of  this  is  given  in  Gear  [8]. 

What  we  intend  to  pursue  here  is  the  case  when  the  step  size  is 
allowed  to  vary.   We  examine  whether  this  will  affect  the  stability  and 
convergence  of  the  method.   We  found  that  with  "mildly"  changing  step  size, 
the  interpolation  method  remains  stable  and  convergent.   A  detailed  proof 
and  discussion  is  given  in  Chapter  3. 

We  substantiate  the  description  of  interpolation  method  with 
an  example . 

Example  1.1 

Three-step  Adams-Bashforth  predictor  and  Adams -Moulton  corrector 
expressed  in  the  form  of  interpolation  method. 


yn,(0) 
hn-iyn,(0) 


n-1     n 
2      yn,(0) 


n-1 


1111 
0  12  3 
0  0  13 
0     0     0     1 


T"yn,(0) 


yn,(m+l) 


hn-l  yn,(m+l) 


n-1     „ 


°'& 


'& 


2      °n,(m+l) 


n-1 


6      yn,(m+l) 


y 


n,(m) 


n-1     n,(m) 


h 


2      yn,(m) 


6      yn,(m) 


(S3 


n-1 


h  v1 

n-2  *n-l 


n-2      „ 


2     "n-1 


n-2 


T"  yn-l 


J 


3/8 
1 

3/U 
1/8 


hn-l   (f(yn,(m)'   tn)   "yn,(m)) 
m  =   0,   ls. . . 


1.2     Variable   Step  Method 

Consider  the   3-step  Adams-Bash forth  predictor  with   constant   step 

size  h 

y     =  y         +  ^  (23  f         -  16  f         +5f     ,)  (1.12) 

n  n-1        12  n-1  n-2  n-3 

One  way  to  derive  (1.12)  is  to  require  that  the  method  be  exact  for  all 

polynomials  of  degree  less  than  or  equal  to  three. 

We  shall  carry  this  technique  over  to  the  variable  step  case  to 

compute  the  coefficients  &  ,  .  §   ,  .  §  n   and  B   -,  in  the  equation 

n,l   n,l   n,2      n,3 


y  =  a    y    +  h  .(§   .  f  -  +  0   0  f    +  B   _  f  ,) 

n    n,l  n-1    n-1  n,l  n-1    n,2  n-2    n,3  n-3 


(1.13) 


2  "3 

Consider  the  polynomials  1,  (t-t  _),  (t  -  t  _)   and  (t  -  t   )  , 

n-3  n-3  n-3 

The   requirement   for  them  to  be  exact   in   (1.13)   yields   a  system  of 


four  linear  equations 

r 


n,l 


=   1 


(l.iiO 


(t     n-t        )   d  +  h  3  +h  fl  +  h  6      ,  -   (t     -t        ) 

n-1         n-3       n,l         n-1     n,l         n-1     n,2         n-1     n,3  n         n-3 

(tn-l   "  W2  *n.l   +  2  hn-l    (Vl   '  W   gn,l 

+  2  h^    (tn_2   -  tn_3)    6n>2    =   (tn   -  tn_3)2 

(t       -  t    .)3  a       +  3  h       (t    .  -  t    _)2  g    . 

n-1  n-3  n,l  n-1        n-1  n-3  n,l 

+  3  h  ft         -  t     „)2  B  =  (t     -  t     J3 

L  n-1       n-2  n-3  n,2  n  n-3 


Upon  solving   (l.lU)  we  have 


r 


a  =  1 

n,l 


3     ,   =  1  + 
n,l 


(1.15) 


hn-l    (2  hn-l  +  6  hn-2   +  3  hn-3) 
6  hn-2    (hn-2   +  hn-3} 


n,2 


h     n(2h          +3h     G+3h      .) 
n-1 n-1 n-2 n-3 

'~  6  h  h 

n-2     n-3 


h  (2  h         +  3  h        ) 

n-1  n-1  n-2 


Pn,3  6  h  (h     0  +  h     ~) 

n-3       n-2         n-3 

This  illustrates   the  derivation  of  variable  step  methods   for  3-step 

Adams -Bashforth  predictor  methods. 

Motivated  by  the  preceding  example,  we  proceed  to   construct  the 

variable   step  method  for  a  general  multistep  method.      For  each  k-step 

method   (1.2)   we   define   an  equation 


E      (d      .    y      .    +  h 


i=0 


n,i     n-i  n-1     n,i     n-i 


f      .  )   =  0 


(1.16) 


where  a      .    and  $   .  are  functions  of  the  step  sizes  h   ,  .  h  ^,..., 
n,i      n,i  n-1   n-2 

h  ,  ,  as  the  corresponding  variable  step  method.  We  require  that 

1.  a.  =  0  =>  d   .  =  0,  3.  =  0  =>§   .  =  0  (1.17) 

i       n,i       i       n,i 

2.  For  some  non-vanishing  a.'s   and  3.'s  we   specify 

the   corresponding  a      . 's   and  3      . 's   such  that 

n,i        n,i 

a   .  =  a.  +  0(h  _),  3   .  =  6.  +  0(h  _)  (l.l8) 

n,i    i      n-1    n,i    i      n-1 

3.  For  a  r-th  order  method  the  total  number  of  prescribed 

&   . 's  and  3   .'s,  vanishing  or  non -vanishing,  should 
n,i       n  ,i 

be  2k  +  l  -  r. 

h.      The  polynomials  1,  (t-t    ),...,  (t-t    )   satisfy 

n  — k.  n  "tK. 

(l.l6)  exactly. 

It  is  worthwhile  to  explain  the  purpose  of  the  above  conditions 

l.-k.      First  let  us  consider  a  linear  difference  operator  L  given  below. 

n 

k 

L  [y(t  )]  =  Z     (d  ,  y(t  .)  +  h    8   .  f(y(t  .),  t  .))     (1.19) 
n    n     .__   n,i    n-i     n-1  n,i      n-i    n-i 

Without  loss  of  generality  one  may  assume  d  .  =  -1.   If  we  substitute 

n,0 

the  Taylor  series  expansion  of  y(t   .)  and  y'(t   .)  =  f(y(t   . ),  t   .) 

n-i  n-i         n-i    n-i 

about  t  ,  for  i  =  1, . . . ,  k  and  factor  out  the  terms  y(t ,  ) ,  y * (t  ,),... 
in  (1.19) ,  we  get 

Ln[^(tn)]  -   C„0  ^n-k'  +  °nl  \-l   *'<W  (1-20) 

where  it:  _  =  1  -  (d  .  +  d  _  +.  .  .+  d  ,  )  (l.2l) 

i  nO        nl    n2       n,k 

J 

\  k  h   .      „   k    k  h 


l!   J-l  hn-l     *l    i-2   J-i  Vl    n'1"1 


LCnq 


k    k  h  . 
E   (  z  -S=J.)<1-- 


{*'l)l   I-l  j=lhn-l      »•" 


10 


It  will  be   shown   in   Chapter  3  that   &      .    and  3      .    are  uniformly  bounded. 

n,i      n,i 

h 

Together  with  the  assumption  that is  bounded,  it  is  clear  that  C 

h  .  n  ,q 

min  '^ 

is  bounded  for  all  n.   When  we  set 


C   =  C  .  =...=  C   =  0  (1.22) 

nO    nl       nr  ' 


we  have 


■kW*."   -   °n;r+l  Cl  ^'r+1)    (Vr>   +   °<Cl>  ^ 

and  this  linear  operator  L  is  said  to  be  of  order  r.   Equation  (1.22) 

can  be  accomplished  in  a  different,  yet  more  convenient  manner  by  condition  k, 

Now  (1.22)  yields  a  system  of  (r+l)  linear  equations  with  (2k+l) 

unknowns  a    , ...,  a    ,  3    ,...,6    .   Dahlquist  [7]  pointed  out  that 
n  ,  J-       n  yK   n , U       n  ,K 

in  the  case  of  constant  step  k-step  methods  the  maximum  order  for  stability 
is 

k+1  if  k  is  odd 

k+2  if  k  is  even 

This  means 

r  <_  k+1  if  k  is  odd 

(1.210 

r  <_  k+2  if  k  is  even 

The  specifications  1.-3.  allow  us  to  arrive  at  a  system  of  (r+l)  non- 
homogeneous  linear  equations  with  r  satisfying  (1.24).   For  such  a  system 

unique  solutions  of  the  unknown  &   . 's  and  3   . 's  exist.   Their  evaluation 

n,i        n,i 

constitutes  the  major  task  of  variable  step  methods.   Therefore,  we  shall 

state  here  an  algorithm  to  determine  all  coefficients  a   .'s  and  3   • 's. 

n,i       n,i 


11 


1.  Impose  conditions  1.-3. 

2.  Use  condition  k.    to  obtain  a  system  of  (r+l)  linear 

equations  with  (r+l)  unknowns  d   . 's  and  3   .'s. 

n,i        n,i 

3.  Solve  for  the  unspecified  d   . 's  and  (3   .'s. 

n,i        n,i 

We  will  prove  that  a  variable  step  method  is  also  stable  for 
"mildly"  varying  step  sizes.   For  Adams  method  the  stability  proof  of 
this  method  requires  much  less  restrictive  conditions  than  interpolation 
methods.   This  advantage  will  be  demonstrated  in  the  numerical  test  in 
Example  2.3.   A  detailed  treatment  of  the  general  proof  of  stability  and 
convergence  will  be  offered  in  Chapter  3. 

For  a  complete  illustration  of  predictor-corrector  methods  in 
the  form  of  variable  step  methods,  consider  . 

Example  1.2 

Three-step  Adams -Bashforth  predictor  and  Adams -Moult on  corrector 
method. 
A-B  Predictor 

y      ,rtl   =  y     ,    +  h_   ,    {&]    f_    ,    +  S.(p>    f_    n  +   &l   f_   J 


where 


n,(0)   "    °n-l  n-1       n,l     n-1  n,2     n-2  n,3     n-3 


/x               h.(2hn+6h_+3h_) 
i(p)   _  -,    +     n-1  n-1 ^2 n-3 

n,l  6  h  (h         +  h        ) 

n-2        n-2  n-3 

/    x  h  (2  h     .    +  3  h  +  3  h        ) 

,IP)   _         n-1  n-1  n-2  n-3 

n,2  '  6  h     _  h     ~ 

n-2     n-3 

r    \       h     .    (2  h     .    +  3  h     _) 
(p)   ..     n-1 n-1 n-2 

n,3  "    Th     0    (h     _   +  h     .) 
n-3     .  n-2         nr3 
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A-M  Corrector 


y  ,   ,  =  y    +  h    (S(c)  f  ,  %  +  6(c)  f    +  3(c)  f    +  B(c)  f   ) 
yn,(m+l)   yn-l    n-1   n,0  n,(ni)    n,l  n-1    n,2  n-2   Pn,3  n-3 


where 


(n)       ,   h  .  (h  _  +  U  h  0  +  2  h  J   h2  1  (h    +  2  h    +  2  h  J 
'Ui  _  1    n-1   n-1      n-2      n-3     n-1   n-1      n-2      n-3 

n,0  ""  2  '      12  h    (h    +  h   _)    '    12  h    h  ,  (h    +  h    ) 
n-2   n-2    n-3  n-2  n-3   n-1    n-2 

h2  ,  (h    +2h  0) 
n-1   n-1      n-2 

12  h  _  (h  "  +  h   )(h    +  h  _  +  h  _) 
n-3   n-2    n-3   n-1    n-2    n-3 

„/  n   n   h  .  (h    +  U  h  Q  +  2  h  _) 
"(c)   1+  n-1   n-1      n-2      n-3 

n,l   2       12  h  "  (h  "  +  h  .) 
'  n-2   n-2    n-3 

t    \  h2.(h,+2h0+2h0) 

-(c)  _    n-1   n-1      n-2      n-3 


n,2  "      12  h    h    (h  .  +  h  ■  ) 

n-2  n-3   n-1    n-2 

(    \                             h2  n  (h    +  2  h   _) 
"(c)  _  n-1   n-1 n-2 

n,3  ==  12  h  TTh     "  +  h   )(h    +  h  "  +  h  Z) 
n-3   n-2    n-3   n-1    n-2    n-3 

f   ,  x  =  f(y   ,  v,  t  )     m  =  0,  1,. . . 
n,(m)     wn,(m)   n 


1.3   Comparison  of  the  Two  Methods 

Both  interpolation  methods  and  variable  step  methods  are  derived 
from  (1.2)  and  they  reduce  to  (1.2)  when  the  step  size  is  held  constant 
throughout  the  numerical  computation.   With  variable  step  methods  one 
computes  y  directly  from  the  k  preceding  points  and  their  derivatives, 
whereas  for  interpolation  methods,  one  computes  y  from  the  Nordsieck  vector 
at  t    .   When  we  represent  the  variable  step  method  in  the  Nordsieck  form, 
we  obtain  a  modified  version  of  the  interpolation  method  with  the  column 
vector  l_  being  a  constant  vector.   Brayton  [3]  has  shown  in  the  stiff  case 
that  these  two  methods  are  not  equivalent  when  the  step  size  varies. 
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It  has  been  definitely  established  in  Chapter  2  that  corresponding 
to  Adams's  scheme  the  variable  step  method  is  more  stable  than  the 
interpolation  method.  For  the  stiffly  stable  methods  the  same  claim  is 
supported  by  numerical  testings  only.   However,  it  will  be  shown  in 
Chapter  3  that  a  more  restrictive  condition,  namely,  that  the  ratio  of 
successive  step  sizes  satisfies 

V^-0^-  (1-25) 

n-1 
is  sufficient  for  either  step  changing  method  to  be  stable  if  the 
corresponding  fixed  step  method  is  stable. 


Ik 


2.      STABILITY  AND  CONVERGENCE   OF  ADAMS   METHODS 

In   this   chapter  we   are   going  to  examine   in  detail  the 
interpolation  method  and  variable  step  method  based  on  Adams -Bash forth 
predictor  and  Adams-Moulton   corrector   (equivalent  names:      Adams   method 
or  Adams-Bashforth-Moulton  method)    methods.      The   convergence  behavior 
of  either  method  will  be  analyzed  and  will  be  exhibited  in  examples. 

2.1     Interpolation  Method 

As   introduced  in   Chapter  1,   the  representation  of  the 
interpolation  method  is   given  by    (l.lO)-(l.ll ) 

Sn,(0)   =  A  °n-l  Vl 

a     i    ^  \  =  a      (    \   +  l    F(a     i    \  )• 
— n,(m+l)        -n,(mj        —       -n,(m) 

With  an  appropriate  l_  this  represents  an  Adams-Bashforth  predictor  and 

Adams-Moulton  corrector.  Let  a(t  )  be  the  correct  value  of  a  at  t  =  t  .   Consi^ 

—  n  —        n 

S   ,nv  =  A  C    a(t    )  (2.1) 

-n,(0J      n-1  —  n-1 

a   ,    .    =   a   ,  x  +  I    F(a  f    v  (2.2) 

-n,(m+l)   -n,(m)   —   -n,(m) 

Then  the  local  truncation  error  d  is  defined  as 

— n 

d  =  a  -  a(t  )  (2.3) 

n   -n   —  n 

Let 

e   /  v  =  a   /  >  -  a,   /  x  (2.U) 

-n,(m)   -n,(mj   -n,(m) 

e  =  a  -  a(t  )  (2.5) 

n   -n      n 

Subtracting  (2.1)  from  (1.10)  and  (2.2)  from  (l.ll)  we  have  (with  the 

notations  given  in  (2.3),  (2.U)  and  (2.5)), 
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e  fns   ■»  A  C  .  e  (2.6) 

-n,(0;      n-1  n-1  ' 


e 

-n 


,  .^e    .  +  Jl  (F(a  ,  J  -  F(a  ,  J)  (2.7) 

,,(m+l)   -n,(m)   -     n,(m)'    vrn,(m)" 

By  the  mean  value  theorem 

F(a  ,  J  -  F(a   ,  J  =  ~  (E   ,  J  e   ,  v  (2.8) 

n,(m)'     -n,(m)    3a  '^(m)  -n,(m) 

where  - —  ( £   /  v )  is  a  row  vector  evaluated  at  a  point  between 
3a_  ■2n,(m)  * 

a   /  x  and  a   /  \.   Substitute  (2.8)  in  (2.7)  to  get 
-n,(m)     -n,(m) 

e  (    J.,\   =    (!  +  £  ¥~   (I   f  J)  e   /  \  (2-9) 

-n,(m+l)        —  3a_  -n,(m)   -n,(m) 

Without  loss  of  generality,  one  may  assume  that  there  are  always  M 

corrector  iterations.   Combine  (2.6)  and  (2.9) 

M-l 
en,(M)    \lQ    <I+i!i  (£„,(!)»   ",-xVl  (2-10) 


Let 


Then 


M_1  3f 

S      =      n      (I   +  *.-g-  (E      ,,J)    A  (2.11) 

n        i=0  3^     -n,(i) 


e      /Mx   =  S     C  e  (2.12) 

-n,(M)  n     n-1  -n-1 


The  local  error  e     given  by    (2.5)    is   rearranged  in  a  computable  form 


as   follows 


%  =  S*  -   a^n)  (2.13) 

=  ^,(M)   "  a(tn} 

=  2^,(M)    "  ^i.(M)   +  ^n,(M)    "  a(tn) 
-n,(M)        -n 


Substituting  (2.12)  into  (2.13) 
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e  =  S   C   .  e  _  +  d 
-n    n  n-1  -n-1   -n 


9F 


For  Adams  methods  —  (£  / . \ )  has  the  following  form 


(2.1U) 


3F 


~  (£   ,,J  =  [h  .  f  ,  -1,  0,...,  0]k 
3a  ^rijCi)      n-1  y  '  ? 

—  — n 


(i) 


=  hn-iV((i)^T-5-i 


(2.15) 


where  £   / . \  is  the  0-th  component  (all  numbering  starts  from  0)  of 
n  ?  \  i ) 

£   / . n  and 

0 


T 
6.  = 

—l 


1+1 


Thus,  (2.11)  becomes 


M-l 


(2.16) 


M-l 


Suppose   f     is  bounded,   i.e.   there  exists   a  constant  L  >   0  such  that 

If  I  <  L 

I     y|     _ 


(2. IT) 


Then  it  is  clear  that  upon  expanding  the  right  side  of  (2.l6)  we  have 


S  =  S  +  h  n  S 
n        n-1  n 


(2.18) 


where 


s  =  (i  -L^M  A 


(2.19) 


and  S  is  a  matrix  whose  elements  are  polynomials  in  h   .  with  bounded 
n  "  n-1 


coefficients.   From  (2.lU)  and  (2.18) 


IT 


e  =(S+h  nS)C  ne    +  d  (2.20) 

-n        n-1  n   n-1  -n-1   -n 

n    n 
*  I         n   (S  C.  ,)  (h,  .  S,  C,  _  e4  .  +  d.) 

+  n\  (s  <W  % 
1=1 

Note  here  the  product  notation 

n 

n 

i=j+l 
is  meant  to  have  the  following  multiplicative   order 

n 

n        (S  C.    J    =  S   C     .    S  C     „...S  C.  (2.21) 

.=.+1  1-1  n-1  n-2  j 

The  matrix  S  which  corresponds  to  a  (k+l) -value  Adams -Bashforth- 
Moulton  methods  for  the  first  order  equation  has  the  following  properties  [8] 
1.   It  has  an  eigenvalue  of  1  and  k  vanishing  eigenvalues. 
Thus, 

Sk+1  -  Sk  =  0  (2.22) 


2.   Suppose  S  =  (s. .).   Then 


3.   There  exists  a  constant  k   >  1  such  that 

s 

llsm|llk        for  all  integer  m  ^  1  (2.2*0 

Let  us   demonstrate   the   above  by  considering  a  U-value  Adams-Bashforth- 
Moulton  method.      From  (2.19) 
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■( 


10      0 

0 

3/8 

0      10 

0 

1 

0      0      1 

0 

3A 

0      0      0 

1 

1/6 

[0,  1,  0,  Q]\ 


M 


1 

1 

1 

1 

0 

1 

2 

3 

0 

0 

1 

3 

0 

0 

0 

1 

1 

0 
0 
0 


3/8     0 

0 

M 

0        0 

0 

3/k     1 

0 

1/6    o 

1 

1 

1 

1 

1 

0 

1 

2 

3 

0 

0 

1 

3 

0 

0 

0 

1 

s2  = 


1 

-3/8 

0      0 

0 

0 

0      0 

0 

-3/1+ 

1     0 

0 

-1/6 

0      1 

1 

5/8 

iA 

0 

0 

0 

0 

-3A 

-1/2 

0 

-1/6 

-1/3 

1 

11/2U 

1/6 

0 

0 

0 

0 

0 

0 

1/1+ 

0 

0 

0 

2/3 

0 

0 

1 

1 

1 

— 

1 

0 

1 

2 

3 

0 

0 

1 

3 

0 

0 

0 

1 

for  any  integer  M  >_  1 


0 

3A 

1/2 
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S3  = 


1  1/2  1/6  0 

0  0    0  0 

0  0    0  0 

0  0    0  0 


s"  = 


1  1/2  1/6  0 

0  0  0  0 

0  0  0  0 

0  0  0  0 

The  fact  that  the  second,  third,  and  fourth  row  of  SJ  vanish  is  highly- 
desirable.   Consider 


2 

U). 

03. 

l 

1 

1 

2 

6 

0 

s3 

c.   = 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Thus, 


3     3       3 

S  C.  S"3  C...  .S°  C  = 
i     J       m 


2 

1  f  f  ° 


0  0  0  0 
0  0  0  0 
0   0   0  0 


(2.25) 
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It  is  clear  that  C  ,  is  an  identity  matrix  when  there  is  no  step  size 

n-1  r 

change  from  t    to  t  .   Suppose  the  step  size  does  not  change  more  than 
every  third  step  and  that  C  i-   I,  then  (2.21)  yields  for  n-j  ^3 

J 


n 

n 

i=j+l 


<S  <W  = 


0  0  0  0 
0  0  0  0 
0   0    0   0 


Furthermore,  suppose  we  set  a  bound  on  the  step  changing  factor  to.  for 


CO,    CO 

1   2 


^ 


0 


all  i.   Then  it  is  clear  that 


:s  Vi* 


i=j+l 

is  uniformly  hounded.   This  is  a  very  useful  result  as  can  be  evidenced 
later  in  the  convergence  proof. 

We  shall  generalize  this  important  result  to  (k+l)-value  Adams 
method  for  the  first  order  equation.   So,  unless  otherwise  specified,  the 
matrix  S  henceforth  used  in  this  chapter  is  based  on  (k+1 ) -value 
Adams  method. 


Lemma  2 . 1 


Let  x  be  a  vector  in  the  real  (k+l) -tuple  space  K 


k+1 


and 


ii  k  T 

x|  |    =1.      Then   S  x  =   vJ±s\   f°r   some   constant  y     which  is  uniformly 


bounded  for  all  such  x. 
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Proof 

Observe  that 

T  T 

T 
lence  L   is   an  eigenvector  of  S  with  eigenvalue  1.      On  the   other  hand, 


apply   (2.22)    on  x 


S(Skx)   =  Skx 


Thus,  S  x  is  also  an  eigenvector  of  S  with  eigenvalue  1.   Since  "both 

T      k  k+1 

6^  and  S  x  are  in  the  same  one-dimensional  subspace  of  R   ,  they  are 

linearly  dependent.   There  exists  a  scalar  constant  y  such  that 

S\  =    yx^  (2.26) 

Now 

UXI«IUX#I  (2.27) 

=    l|Skxl| 


iiiskn  iiiii 

k 

—         s 

So  we  see  that    |  \j    \    is  uniformly  bounded  for  all   |  |x|      =  1. 

Q.E.D. 


Lemma  2 . 2 

h 

Let  W  >  0  be  a  constant  such  that <  W.   (Note  that  it  is 

min 
h 

natural  to  assume >  1,  •.  W  >  l) .   Suppose  the  step  size  does  not 

min 

change  more   than  every  k  consecutive   steps.      Then 


N 

n      (sc.  J 

1=J+1 

is  uniformly  bounded  for  all  j  =  0,  1,...,  N-l. 
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Proof 

First  we  observe  that 

for  an  arbitrary  step  changing  matrix  C. .  Hence,  as  long  as 

N  k 

n        (S  C        )    contains   a  factor   S    ,    (2.26)    and   (2. 27)   lead  to 

1=3+1  X 

N  N 

||     n      (s  c      )|  |  =     sup      ||     n      (s  c      )  x|  | 

i=j+l  11x11=1       i=J+l  X 

<  k 
—    s 


Let  the  J  be  the  index  set 


N 


{The  values  of  j  such  that   n   (S  C._  )  fails 

i=j+l     X 

k 

to  contain  a  factor  S  } 

From  the  hypothesis,  J  contains  less  than  2k-l  elements.   Thus,  we  see  that 

N 
Max  ||   n   (S  C.  _)||  i  (k  ^fk~2 

v   2k— 2 
Let  k     =    (k     w    )  .      Summing  up  the  above  two   cases,  we  have 

J_  s 

N 
II      n       (S  C.      )||    <.k  if.  3  (2.29) 

i=0+l 


The  bound  k  depends  neither  on  N  nor  on  j . 


Q.E.D. 


Since  we  are  dealing  with  methods  of  at  least  second  order,  the 

local  truncation  error  d.  satisfies 

i 
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ijl  <  h._x  e   ¥i  (2.30) 


where  e  is  0(h.  _).   It  is  also  clear  that  there  exists 

l-l 

a  constant  k     >   0   such  that 

||S     C        ||    <  k  ?j  (2.31) 

Having  established  bounds  for  various  terms  in  (2.20),  we 
are  in  a  good  position  to  present  the  convergence  theorem. 

Theorem  2.1 

The  interpolation  method  based  on  (k+l) -value  Adams-Bashforth- 
Moulton  method  for  the  first  order  equation  (l.l)  is  convergent  under 
step  size  changing  conditions  provided  the  change  is  not  made  more 
frequently  than  every  k  steps. 


Proof 


At  time  t„  the  error  e,T  satisfies 


N    N 


Let 


Then 


e,T  =  Z    n    (S  C.  .  )(h.  _  S.  C.  .  e.   +  d.) 
-N   J=1  i=(.+1     i-l   J-l   J   J-l  -J-l   -J 

N 

+  n  (S  C.  .)  e. 

.  .     i-l  —o 
1=1 


E 


x  =   |e  I   +  ^-  (2.32) 

n   ' ■  n ■ '   k2 


Sllil-IA  ."     (s  ci-i)(hj-i  §j  ViVi  +  ij)|l  +  H"    (sci-i'Sol 

J      _L     -L      J      J- 

N 
-^    kl  k2  hj-l  Xj-1  +  kl    l|eo'l  (2-33> 


2k 


by  applying  (2.29),  (2.30)  and  (2.31).   Consequently, 

*i  *  ki  k2  hj-i  xj-i +  S  N-oll  +  ir-  <2-3u> 

We  shall  show  "by  mathematical  induction  that 

k  k   (t„-t  ) 
^i^ll-oll*^)-1  2   "  °  (2-35) 

1.  When  N  =  0,  consider  (2.32) 

<^II«0II  ♦  £ 

since  k  >_  1. 

2.  Assume  (2.35)  holds  for  N  £n 
From  (2.3*0 

n+1 

x  j.1     :     z     k     k.   h.        x.   .    +L    i|eji    +f 
n+1  —  1     2     j-1     j-1         1    '  '    0'  '        k2 

n+1  p         k.   k      (t.      -t_) 

j=l  d 

+  k,     lie'   II   +  .— 

1    M    0M       kg 

n+1  k     k0    (t.      -t    ) 

=   (kl||e0||    ♦  £)<!  ♦     ^  ^  ^  Vl  .     ^     2       ^     °) 
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The  last  factor  is 

n+1  kl  k2  (ti-l-tO) 

1  +  E  k  Is.     h    e        J--L  0 

j=l       J 

n+1     k     k     (t        -t    )  n+1     k     k     (t        -t    ) 

=  1+Ee12        Jl°(l+kk     h        )   -     Z      e  1     2        J_1     ° 

j=l  J  j=l 

n+1     k     k      (t.-t    )        n+1     k     k      (t        -t    ) 
1l+Le12        J°-Ze12        J_1° 

kl  k2  (WV 
=  e 

Therefore, 

k  k   (t  -tj 
,.   it   ||    e  n    1  2   n   CT 
x  ._  <  (k    e0   +  ■ — J  e 
n+1  —   l  '  '  0M    k? 

The  inequality  (2.35)  is  verified.   It  immediately  follows  that 

Me  II  <  (k   Me  II  +-M  e1"1^  ^"^  -  -£-  (2'36) 

l|eNM  1  Ux  ||e0||  +  k^j  e  ^ 

e  is  of  order  0(h    ).   Assuming  the  use  of  Runge-Kutta 
max 

method  or  Euler's  method  for  the  initial  points  needed  to  start  the 

multivalue  method,  we  must  have   leM  I  ■+  0  as  h    -*  0 .   Hence,  as 

11    0 ' '  max 

h  ->  0,       le^l  |    ■*■  0   for  all  N  such  that   t.T  e    [0,b].      This   means   the 

max  ' '    N1 '  N 


computed  solution  converges   to  the  exact   solution  as  h  ■+  0 

max 


Q.E.D. 


In  the  preceding  proof  the  step  size  change  is  restricted 

to  every  k  steps  in  order  that  n  (S  C.   )  based  on  a  general  (k+l) -value 

i 
Adams  method  would  be  uniformly  bounded.  We  like  to  point  out  that  this 

restriction  could  be  removed  in  the  case  of  the  3-value  Adams  method 
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I  = 


5/12 

1 


1/2 


S   = 


10   0 

5/12 

0  10 

- 

1 

0   0  1 

1/2 

jo,  i,  oJ 


— 

1 

1 

1 

0 

1 

2 

0 

0 

1 

Consequently, 


1     7/12     1/6 
0        0  0 

0      -1/2        0 


for   any  integer  M  >   1 


1 

7uJ 

12 

2 

ID...      03. 

_   J+1    J 

12 

1  2 

0 

0 

0 

0 

0 

0 

N 

n      (s  c      )  = 

i=j+l 


Since  o> .  <  W  ¥■  j,  the  above  matrix  is  uniformly  bounded.   Therefore,  we 
J 

can  obtain  convergence  without  being  restricted  on  the  frequency  of  step 
size  change. 

However,  the  situation  is  different  for  higher  value  Adams  method, 
If  the  step  size  change  is  subject  to  gross  changes,  then  the  error  may 
grow  rapidly  throughout  the  numerical  computation.   Consider 
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Example  2.1 

The  error  growth  of  the  U-value  Adams  method  with, 
alternating  step  size. 

In  this  example  we  are  mainly  concerned  with  the  behavior  of 

N 

n   (S   C        ).      It  is  easy  to  see   that        II        (S   C        )    can  be   converted  to 
i  1_  i-j+1 

the  matrix 

1       X  X  X 

0     0  0  0 

_    .    .    TO    .N-l   _   N-l  M    ...    ,      .    _    ,    _      N-l   0  N-l 

.  /    _vN-j,lvN-j  2  ,  s  f       xN-j+l,3wlxN-j-l  2     _      ,  » 

o    x    (-1)    °(— )        n  w.     n    (l-w  )  (-D         (r)(p)         to    n  u.     n    (i-a>.) 

1=0      i=j+l  °i=j      i=j+l 

o    x    (.dN-^^^-mV^  Yu-uO     (-i)^+1(^u.Y^Y(i-u.) 

J     *  i=j  1i=j+l         x  Ji=j     i=j+l       x 

where  x  denotes  other  non-zero  elements  in  the  matrix.   The  eigenvalues 
of  this  matrix  are  determined  to  be 

„  .  N-l  £    (l-w.  ) 

i,  o,  o,  (-i)N"J  n   x  2  x 

i=j 

Of  the  above  eigenvalues,  only  the  last  one  is  affected  by  step  size  change. 

Suppose  for  any  two  consecutive  steps  the  product  of  their  step  size 

changing  factor  is  unity.   Then  for  some  oj  >  1 

-1 
w2j  "  w2j+l  ~  w 

When  N-j  is  even 
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2 

H    ,    N-l    w-    (.1-0)   )  N-l     1-w 

(-DN"J     n       X    §    *   1  -  I  n     (-ylll  (2.37) 

4o) 


when  N-j    is   odd 


_    .   N-l   ok     (1-oa.)  ,_       ,    N-2   o)2    (l-o).) 

|(-DN-J     n    JL_^.pIl=al    n      i         1  1  (2.38) 

i*J  i=J 

.   u,(a,-l)    r(o)-l)21    ^i 

-  ~2     [-^r~]    2 
(  -l)2 

It  is   clear  that  — f is  an  increasing  function  of  o)  for  o)  >   1.      Thus, 

4o) 

(o)-l)2 
it  is  of  interest  to  find  the  value  of  o,  >  1  such  that  — j =  1,  i.e. 

40) 

(o)-l)      =   4o).      The  roots   of  this   equation   are  3  +  2/2.      This   implies   that 

(I)2 

VU),  >   1  when  o)  >    3  +  2/2  (2.39) 

4oi 

Hence,  for  a  relatively  large  u)  the  magnitude  of  the  eigenvalue 

n-i  a,2  d-o).)     n 

(-1)  II of       n        (S   C        )    grows   rapidly  as   N-j  becomes   larger. 

i=j  i=j+l 

Consequently,  the  error  er  grows  rapidly  according  to  (2.20). 


Example  2.2 

Numerical  testings  of  the  4-value  Adams  method  with  condition 
(a)  or  (b); 

(a)  step  size  being  alternated 

(b)  step  size  being  alternated  once  every  three  steps 
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Numerical  testing  "with  conditions    (a)    and  (b)    are  performed 
on  the  following  two  problems  using  the  IBM  360/75. 

1.  y'  =  -y       y(o)  =  1  (2.U0) 

2.  y1  =  Ut3   y(0)  =  0  (2.Ul) 
The  initial  step  size  and  step  size  changing  factor  are  taken  to  be 

0.05  and  10  (or  O.l),  respectively.  We  are  interested  in  comparing 

computed  solution  and  exact  solution  throughout  the  time  interval  [0,5] • 

We  summarize  such  results  in  Table  2.1,  Table  2.2,  Table  2.3  and 

Table  2.4.   (A  representative  program  (for  Table  2.2)  is  given  in  Appendix  A.) 


TIME 


COMPUTED 
SOLUTION 


EXACT 
SOLUTION 


0. 72? 5 

C.  5488 
0.4168 
0.3166 
0.240b 
0.  18  26 
0.  1387 
0. 10  53 
0.8005 
0.oT81 
0.4618 


C.335 
0.362 
C.390 
0  .4  17 
0.445 
0.472 
0.50C 


0. 3508 
C. 2664 
C.20  24 
0.1537 
0.1167 
0.8870 
0. 6  73  7 


2  73536D 
1  163610 
6201970 

3  676940 
084632D 

8  3  52410 
613L220 
9 92 2 46 D 
8312790- 
0062630- 

9  6  2333  D- 
*35410~0- 
9097340- 
191  145  0- 
5191660- 
8566970- 
7139100- 
94b  9990- 


00 

00 

oc 
66 

Tie* 

00 
00" 
01 
01 
_0_1_ 

n 

01 

Jl 

01 
01 

02_ 
0  2 


Table  2.1.     Integration  of  y'    =  -y,  y(0)   =   1  using  U-Value  Adams 


Interpolation  Meithodwwith  Alternating  Step  Sizes   of   .05  and   .005 
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TIME 

0.28000D 
0.6CCCCD 


COMPUTED 
SOLUTION 


EXACT 
SOLUTION 


00  0. 75578385720    00 

00  C. 548811 «8770    00 


0. 


75578 
5AB_ai 


374150  CC 
16361C  QO 


O.87500D  00     0.41686234420  CC     C. 4 1 6 8 e 2C 1 c 7D  CC 

C.11050C  CL 0..J31211186  5D  DJ) 0.  33  121C  ee22C  CCu 

0.142500  01     C.24C5C87576C  CC     0. 24C5C846320  CO 

._G.17000C  01 O*.l£26_8J?£0eJ.D_GC C.  1826  82524  ID  _GC 

0.193C0O  01     0.  1451484429C  00     0 .  1 4 5 148  1985C  00 

0.225CQD  ci    QJL££52S2£2*t2n   c.c          c.  i o 57qq?74<sn  on 


0.252500 
0.2755  0C 
0.3C75CO 
0..33500C 
0.358C00 

c.3qccoo 


CI 
_Q1_ 

CI 

OL 

01 

iLL 


C. 80058497660-01 
_CL.6J6C5i7  54  8C-0L 

C. 46189755630-01 
JL^3_50  844615  8C-01_ 


C. 27875787870-01 
c.?r?4i<;pp7cn-ni 


c. 

0. 
C 
Q. 
0. 
_C^ 


eOC58 
63609 
4618c, 
3508  A 
27875 
?  C  '  4  1 


212790-C1 

G1967C-CL. 

628380-01 

2541CD-CL 

69826C-01 

^1  145P-C1 


0.417500    01 
L_Cl*A4Cl5J)D_C_ 
0.47250C    01 
3£Q0H_Q1 


O.153752507OD-01  0.15275 

.  Z2±tl  13  C  4  C  -  CJ 0*  12  21 6 

0.8e7C751928D-02     C.8e707 
L^7  3797_a032  C-Q2 0.67379 


1S166C-C1 
1C641E-CL- 
1291GD-C2 
6  95L9  £r_Q2_ 


Table  2.2.      Integration  of  y'   =  -y,  y(0)   =  1  using  H-Value  Adams 
Interpolation  Method  with  Three  Steps  of   .005  followed  by  Three 
Steps   of   .05 


TIME 

0 .325000 
0.60CC0O 
0.875C00 
0.115C00 
0. 1.42500 
O.17C000 
0.197500 
C. 225000 
0.2525^0 

o.?8or^o 

307^00 
3  3  5  C  0  0 

367S0D 

39rcoo 

417  5  0  0 
0.445000 
0.  472^0 
0.500000 


00 

00 

^0 

CI 

CI. 

01 

01 

01 

01 

01 

01 

CI 

01 

01 

01 

01 

oi 

01 


0. 
0. 
0. 
0. 
■0. 
C. 

~    • 

0. 

■c. 


•o. 
c. 

.n 

v>  • 

0. 

,  n, 

r. 


COMPUTED 
SOLUTION 

11 15500  7440-01 
129  91923140  CC 
575R7097^7P  00 
2  1009461440  ~n\ 
7f5  918f0130__0.1. 
41636942660  03 
138  7_8P2996_D.  05 
4  7 3 C 993 5 530  C6 
161PPA1726P  C8 
548cn3i34on  nq 
1867*^7 6 ?CO  11 
6  3  5  9  5  P  3  8  3  3  0  1  ? 
2\(  54^16^80  14 
73736013720  IS 
2  5107  57^4  10  17 
85452876730  '1R 
2f">1  l.C  8  62  2  00.  ,.?_0 
99124  31664D    21 


_o. 
o. 

0. 

0. 

.p.. 

•  '  . 

.0.. 
0. 


0. 
C. 
C. 

0. 

_o_. 

0. 


EXACT 
SOLUTION 

1 1.1.5  6  64  0620-01 

1296C0C0O0O  00 
5861 81 64  06  0  ^0 
17490C625C0  01 
4,1234378919  .01 

8  3 521 rroono   31 

15714? 7S 390  02 
2562890625O  02 
40648594140    02 

61  46c;6'\0C0O  H2 
81408844140  0? 
12594^06  ?o  0  3 
177676^7540  03 
23134410000  03 
3^787663790  03 
3^213900670  03 
4Aa|4  3  35716P__07_ 

62  50C  000000    03" 


Table  2.3.      Integration  of  y'   =  H   ,  y(.0)   =  0  using  U-Value  Adams 


Interpolation  Method  with  Alternating  Step  Sizes  of   .05  and   .005 
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COMPUTED 

EXACT 

TIME 

SOLUTION 

SOLUTION 

0.28000D 

00 

0.61506441670- 

•02 

0. 

614656CCCCD- 

■C2 

n.ftrrron 

()() 

n.  i  ?9M7?3i9n 

nc 

Ot 

i  pqisrnnnnnn 

no 

0.8750C0 

CO 

C.5662C21C1C0 

oc 

c. 

5e6l6164C60 

00 

£..11050C_CL. 

0^149  C926.55.4CL 

_CL_ 

C1.45C9C2C51D 

XL 

0. 1425C0 

CI 

C.412347C543C 

01 

0. 

,41234378910 

01 

_Q.O70aQD_Cl 

0.P3  5214C8.7  9Q_.C1_ 

C. 

-£3521CC0G.Q0. 

OL 

0. 193C0C 

01 

. 0.13874924930 

02 

0. 

136748 8CC10 

C2 

n.  ??«;nnp 

m 

n.?^f  ?RqR<n?n 

n? 

Or 

?Sfppqn^?^r 

0? 

0.2525C0 

CL 

C. 406  4  86  5  5440 

02 

c. 

4C€485C4140 

02, 

_0.2755QC 

.01 

IU5X608544  84Q. 

J02__ 

0. 

5760E4?9500_C 

0.3C75CO 

CI 

0.  e<54Cfi91763D 

02 

c. 

69408844140 

02 

_Q.335aaa 

01 

0^125944  58800 

_03 

0* 

1259445C62D 
16426C1090C 

-C3- 

0.3580CD 

CI 

0.1642601947D 

03 

0. 

03 

n. ^qnrnn 

ri 

n.  ?^i  iaai  ^cr 

r.-K 

n. 

?-*1  ^Ainrrr 

P3 

0.417500 

01 

0.30382679CC0 

C3 

c. 

30382668790 

C3 

Q.44C5flrL 

.01 

CI 

r,i765ifc?q?4C   n 

0^ 

.37651618620 

-03- 

0.472500 

0.49P433646C0 

C3 

0. 

4984335316D 

03 

0  .  50DQQXL 

-01 

0.62500012260 

03 

625CC0C-QC QD    03 

Table  2.k.      Integration  of  y'  =  Ut-3,  y(0)  =  0  using  U-Value  Adams 
Interpolation  Method  with  Three  Steps  of  .005  followed  by  Three 
Steps  of  .05 

It  is  clearly  demonstrated  in  the  above  example  that  condition 
(b)  is  a  stabilizing  condition.   Naturally  our  next  objective  is  to 
generalize  Theorem  2.1  to  cases  other  than  Adams  methods.   In  view  of  the 
fact  that  the  crucial  condition  (2.22)  is  lost  in  so  doing,  some  modifi- 
cations on  the  hypothesis  of  the  theorem  are  necessary.   This  attempt 
eventually  leads  to  a  much  bigger  gain.   Indeed,  we  are  able  to  prove  the 
convergence  of  interpolation  methods  based  on  .general  multistep  methods, 
of  which  Adams  methods  are  only  particular  cases.   We  shall  offer  the  proof 
in  Chapter  3.   Basically,  it  calls  for  mild  step  size  change  consistent  with 
the  implication  of  the  hypothesis  in  Theorem  2.1. 
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2.2  Variable  Step  Method 

Consider  the  fixed  k-step  Adams-Moulton  method 

..+  3,.  t  .  ) 


(2.U2) 


yn  =  yn-l  +  h(3Ofn+3lfn-l+-"+ek  Vk 

The  corresponding  variable  step  method  is 

u  v,    (a  f  +  6    f    +-.-+  i  ,  f  J      (2.U3) 

yn  =  V!  yn.x  +  hn_l(Bn,0  fn   3n,l  n-1        n,k  n-k 

according  to  part  1  of  the  algorithm  stated  in  Chapter  1.   Now  part  2 
leads  to  a  (k+2)  X  (k+2)  system  of  linear  equations 


"  n,l 

(Vl-U)an,l  +  VlVo+-+Vl^k 


(t. 


\2   *  o  v,  (  +      -  t        )    B  +...+  2  h     ,    (t 


=  1 

=  <*n 

-   *      J 

n-k 

!  Vi 

(tn-k+l  • 

■  tn-k)    Sn,k- 

=  (t    -  t     v 

n         n-k 


(2.UU) 


k-1  "  V/+1  V  +    (M)   Vl   Un  "  V/   *n,0  +  '"  + 


-  t     , .) 


k  s 


=   (t     -  t     .  ) 


k+1 


v. 


(k+D   hn_!   ^n_k+i  "  Vkj      Pn,k-1  "    "n  "     n-k 


This   is   easily  seen  to  be 

C  A 


5n,0    +  ^n,l+...+     n,k 

..+   2    (t 


=  y. 


2 ■«*„  -  Vk'  VO  +  2  (Vl  "  Vk>  Bn,l  +' 


n-k+1       Vk'   Sn,k-1 


( 


=  u, 


(2.U5) 


<k+1>    <*n  "  V/  §n,0  +   (k+l)    Un-1  "  Vk}      *n,l 

^  +  ...+   (k+1)    (tn_k+1  -  tn_k)k  3n>k_1  =  Vl 


33 


where 


y.  = 

1 


(t   -  t  J1  -  (t    -  t  v)i 
n    n-k      n-1    n-k 


n-1 


Let  M.  be  the  matrix 

A 


2  (t   -  t   ,  ) 
n    n-k 


2  (tn-l-  W   '• 


2  (t  .  ._  -  t  ,  )    0 

n-k+1    n-k 


<**><*» -W11  (k+1)  (tn-i  -  Vk)k   •••   |w)(WrVk»k  ° 


and  M.  (i  =  0,  1, . . . ,  k)  be  the  matrix  obtained  by  replacing  i-th  column 

of  MA  with  the  column  vector 
A 


*k+l 


We  claim  that  M  is  non-singular  since 

A 


k-1 

|M  I  =  (-l)^2  (k+i)!  n  (t   -  t  .  ) 

A  ._0  ~  n-x    n-k 


(t   -  t  .  )     ...   (t  ,  ..  -  t  ,  ) 
n    n-k  n-k+1    n-k 


(t  -  t  v)k_1   ...   (t  _._  -  t  v)k_1 
n    n-k  n-k+1    n-k 


^ 


Vandermonde  Determinant 


3>* 


Therefore,  it  is  clear  that  6   .  can  be  uniquely  solved  by  the  formula 

n ,  l 


0 


M. 

l 


n,i    |M, 


i  =  0,  1,...,  k 


(2.U6) 


To  make  (2.U6)  accessible  for  further  analysis  we  need  to 
express  M  in  an  alternate  form,  i.e. 


M. 


k 

2  E   h 

i=l 


n-i 


2  I      h   . 

i=2  n"X 


k  k 

(k+D(  z   h   .r   (k+i)(  e   h  .  r 

i-1  n_1  i=2  n-X 


2  h 


n-k 


(k+l)  h*  ,      0 
n-k 


Then 


_  ,  (i+2+...+k) 


MAi  =1^:k 


i 

k  h 


2  E 


n-i 


1 
k  h 


2  E 


n-i 


i=l  n-k 


h 


i=2  n-k 


i=l  n-k 


h  , 
i=2   n-k 


1    1 


2    0 


k  h  k  h  . 

(k+l)(  E  lT^r  (k+l)(  E  -^r   ...   (k+l)  0 


Now  let 


6i=   h. 


h.  -  h. 
i    l-l 


i-1 


(2.1*7) 


Equivalent ly. 


1  +  6.  = 


i   h 


h. 

l 

i-1 


(2.1*8) 


(2.51) 
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Note  that  6.  may  posses  negative  values.   For  i  <  k 

n-i  _  n-i       n-(i+l)        n-(k-l)    _  /n    .    N    ,   ,  . 
h h "5 "•  ~h n  (i  +  <s  .)    (2.1+9) 

Vk   V(i+l)   hn-(i+2)         hn-k     3-1       n"J 
Thus,  |maI  is  a  (k+l)(k+l)  determinant  whose  (p,q)  entry  is 

1  for  p  =  1 

=  1  °  _         k-1  k-1  for  p  >  1,  q  -  k+1 

P^l     (1+   E    n   (1  +  6n  i^P   f°r  P  >    1.  ^k+l 

I-q  J=i        J  (2.50) 

The  determinant  I M. I  can  be  converted  into  a  similar  form. 
1  i1 

Hence,  it  is  clear  from  (2.k6)   that 

g    _Pi(4n-l'  6n-2"--'6n-k+l) 
n'i   Wl>  6n-2'"->  Sn-k+l> 

where  P.  and  Q.  are  polynomials  in  6   -,,6   ^,...,6  , 

l      l  n-1   n-2       n-k+1 

6    =  6    =...=6      =  0  ¥  n 
n-1    n-2       n-k+1 

the  method  reduces  to  the  fixed  step  method.  We  have 

P.(0,  0,...,  0) 

1  _  Q 

Q1(0,  0,...,  0)  "   i 

This  means  $   .  is  a  continuous  function  of  6   ......  6   ,  ,.  in  a 

n,i  n-1  n-k+1 

neighborhood  of   (0,   0,...,   0).      Thus,   for  some   constants   d  and  D  there 
exists   a  compact  region 

U  =  {(6n-1 Ul'I^Vl Sn-k+llD    *} 

of  the  real  (k-l)-tuple  space  on  which  3   •  is  continuous.   Moreover 

n,i 

I g   . I  is  uniformly  bounded  in  U .  We  shall  state  the  above  results  in 
n,i' 

a  lemma. 


When 


(2.52) 
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Lemma  2 . 3 

The  coefficients  &   .  *s  and  3   . 's  in  (2.1+3)  are  uniformly 
n  ,i        n  ,i 

bounded  provided  for  some  constants  d  and  D  the  step  size  change 
satisfies 

d^„.i WiD  *"  (2-53) 

Now  we  proceed  to  prove  stability  of  the  method.   The  next 
lemma  is  essential  for  that  prupose. 

Lemma  2 . k 

Consider  the  linear  difference  equation 


k 

&_  =  z„  i  +  h„  i  (  z  3_  a    z„  •  +  O  (2.5*0 

i=0 


Jm    m-1    m-1  _._n     m,i  m-i    m' 


Let  B,    q,    A,    z  be  non-negative   constants   such  that 


k 

E      |3      . I    <  B     ¥  m  (2.55) 

1=0       m'1     " 

h      t    3      n   <.<!  <   1     ^m  (2.56) 

m-1     m,0 

|A    |     :   A     ¥  m  (2.59) 

1    m'    -r 

|z.|    <_  z,   0^i£k-l  (2.58) 


Then 

B 


z    I    <  e1"*       m         k~1     (z  +  -^-(t     -  t        ))  (2.59) 

m1    —  1-q        m         k-1'' 
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Proof 

Rearrange    (2.5*0 

z     -h     n    3      _    z     =z         -  h     _    B      nz         +  h     .    3Z   0    z    n 

m         m-1     m,0     m         m-1  m-1     m,0     m-1  m-1     m,0     m-1 

k 

+h     .(    Z      3      .    z      .    +  A    ) 

m-1    .  m,i      m-i  m 


h 

m-1 
.*.   z     =   z 
m  m 


i    +  i VT   "fi ((3      n  +   3      .)    z      _    +      E      3      .    z      .)         (2.60) 

-1       1  -  h     .    3      n  m,0  m,l        m-1        .    _     m,i     m-i 

m-1     m,0  i=2        ' 


Let 


h 

+  2=1 x 

1  -  h     ,6     n     m 
m-1     m,0 


xQ=    |*|  (2.61) 


Then 


z      <      z 


x .    =  max  {  I  z .  I  ,   x .    .,  }        i  =   1 ,   2  , 
l  '    i '         l-l 


'((le.nl  +  K  ,|)U_  J   +    I     \B     Jl*    J) 


m1  —   '    m-11         '  1  -  h      .    3      n  m,0 '         '    m.i '  '    m-11         .    ~    '    m.i '  '    m-i 

m-1     m,0  i=2 

h      i 
m-1  ,    ,       , 


1  -  h      ,3      n  m 

m-1     m,0 


h     _  k  h     _ 

I m-1 I  /   _  I  0         I  \  ,  I m-1 I     I  ,     I 

-  x     i    +     ~i Z o (    l  3      •    )   x  + A 

—  m-1        '  1  -  h     n3      A'    .    A  '    m,i'        m-1      '  1   -  h      n    3      -.      '    m1 

m-1     m,0     i=0  m-1     m,0 

<    (1  +-f^— )    x  +-2izl   A  (2.62) 

—  1-q  m-1  1-q 


Certainly 


h      i    B  h      i 

x  <    (1  +  -^— )    x1+^A  (2   63) 

m-1  —  1-q  m-1  1-q  v^.<~o; 
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Now   (2.61),   2.62)    and  (2.63)    lead  to 


h     .    B  h     . 

xm  1  d  +  —=r~r, — )x™   i    +  TT  A  (2.64) 

m  —  1-q  m-1  1-q  ' 

h  B/l-q  h      . 

m—  1.  m— 1 

—  m-1  1-q 

Bh     . 

m~l  Al 

=  e  x     .    +  Ah     , 

m-1  m-1 

R  A 

where  B  = ,   A  =  ■= .     We   claim  that 

1-q'     1-q 

xmle        ^   (Z  +  A(tm-Vl))  (2.65) 

This  is  verified  by  mathematical  induction 

1.  m  =  0;  obvious 

2.  Assume  (2.65)  holds  for  m  <_n.  From  (2.64)  and  (2.65) 

Bh 

x  _,_,  :  e  n  x  +  Ah 
n+1  —      n     n 

Bh   B(t  -  t,   ) 

:  e  n  e   n    *_1   (z  +  A(t  -  t  - ) )  +  Ah 

—  n    k-l       n 

=  e   n      k"1  (z  +  A(tn  -  tk-1})  +  Ahn 

B(t  .-  -  t.  .)  B(t  ._  -  t.  _) 

n+1    k-l   /   ,  ■*,.     ,    *»       n+1    k-l  ~ 
<e  (z+A(t-t,_))+e  Ah 

—  n    k-l  n 

B(t  ..  -  t,  .  ) 
<e   »*    W  (.  +  A(tn+1  -  Vl)) 

Thus,  the  induction  is  completed.   (2.6l)  and  (2.65)  yield 

B(t  -  t,  ,  ) 

z   <  e   m    k_1  (z  +  A(t  -  t  J) 

m  —  m    k-l 


B   (t  -  t,  J 


=  e 


1-a'm   k-i(s  +  JL(  t  -tvl)) 

1-q   m    k-l 

Q.E.D, 
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We   say  a  variable  k-step  method  is   stable   if  the  sequence 
{y  },  n  =  0,   1,   2,...   is  uniformly  bounded.     We  are  ready  to   show  the 
stability  of  the  method. 


Theorem  2.2 

A  variable  k-step  Adams-Moulton  method  is   stable  provided 
condition   (2.53)   is  satisfied. 

Proof 

Let  us  write    (2.^3)   in  the  form 

y     =  y         +  h  (8      a    (f(y       t    )   -   f(0,  t    ))   +...+  (2.66) 

n  n-1  n-1       n,0  n,n  n 

8     .     (f(y     .  ,   t        )    -  f(0,   t        )))   +  h     .    A 
n,k  n-k       n-k  n-k  n-1     n 

=  y         +  h  (5         f'  (n        )  y     +...+  8         f    (n     .  )  y     ,  )   +  h     .    X 

n-1  n-1       n,0     y      n,0       n  n,k     y      n,k       n-k  n-1     n 

where   X     =  8  f(0,   t    )   +...+  0  f(0,   t     .  )  (2.67) 

n  n,0  n  n)k  n-k 

and  n   .  is  an  intermediate  point  between  0  and  y   . .  By  hypothesis 
n,i  n-i 

in  (l.l)  ,  r^  T.i|f(0,t)|  is  bounded.   Let  B  be  the  uniform  bound  for 
te  [0,b] '    '   ' 

8   .  ,  i  =  0,  1, Let  A  =  (k+l)B  ,  ^\i\  f(0,t)|  .   Then 

nsl  t  LUjDj 


|X   <  A  ¥  n 
1  n1  — 

It  is  clear  that  with  Lipschitz  condition  one  can  choose  appropriate 

h    to  get  (2.56)  for  some  0  <  q.  <  1.   Suppose  we  denote  8   •  =  8   •  f  (n   .) 
n-i  n,i    11,1  y  n>i 

Then  (2.55)  is  fulfilled  if  we  put  B  =  (k+l)  B  L.   The  remaining  condition 
(2.58)  in  Lemma  2 .k   is  trivial,  just  take 


Uo 


v    -  max  I 

Y  ■  v 

0<i<k-l    |yi 


Apply  Lemma,  2 .  1+  to   get 


|y  I   <  ex-q      n        k_1  (T*rfr(t    -  t.    .))  (2.68) 

1    n1    —  1-q       n  k-1 

Bb 


:  el-cL   (Y  +^b_}     ¥n 
-  1-q 

Q.E.D. 

Convergence  follows  immediately  from  stability. 

Theorem  2 .3 

A  variable  k-step  Adams -Moulton  method  with  step  size  change 
satisfying  (2.53)  is  convergent. 

Proof 

Consider  the  associated  difference   operator  L     for   (2.1+3) 
LQ(y(tn)]  =  y(tn)  -  y(Vl)  -  Vl  (BnjQ  f(y(tn),  tn)  ♦...*  (2.69) 

~  ,k+2      (k+2)    f      v  ,.  » 

=  C     .  ._  h     _    y  (t    )        t     e    (t     .  ,  t   ) 

n,k+2     n-1  n  »        n  n-k       n 

since  (2.1+3)  is  of  (k+l)-th  order.   Define  e  =  y  -  y(t  ).   Combine 
(2.1+3)  and  (2.69)  to  get 


1(1 


en  "  en-l  "  Vl    (in,0   (f(V   V   "   f<r<V'   *n»   +  "'  +  (2-70) 

g„,k   (f(y„-k>   W"   f(y(Vk''   WJ) 

+  h  <-C     ,+?hkVk+2>    (T    >) 

n-1  n,k+2     n-1  n 

"  Vl  "  hn01   (§n,0   fy<?n,0»    Sn  +" "  "+   V   V^n,*'    'n-*' 

+  h    .  (-c    k+9^Uik+2)  <*  )> 

n-1  n,k+2     n-1  n 

where  5   .  is  an  intermediate  point  between  y  and  y(t  ).   Assume  y      is 
n,i  n        n 

bounded  in  [0,b].   Then  it  is  obvious  that  Lemma  2.h   is  applicable  here 
when  we  set 

6   .  =  3   .  t  J  C   .  ) 

n,i    n,i   y  n,i 

,   _  a      v,k+l   (k+2),   , 
An  "  -Cn,k+2  hn-iy     (xn} 

As  analogous  to  the  proof  of  Theorem  2.2,  given  some  0  <  q  <  1,  one  can 

choose  some  suitable  h  .  to  satisfy  condition  (2.56).   For  the  other 

n-1 

conditions   in   Lemma  2.U,  we  need  to  set 

~  ,k+l       max      I    (k+2)    ,    .  , 
X*      A  =   Chmaxt£[0,b]ly  0(t)l 

where   C   is  the  uniform  bound  for  C     ,  ._ 

n,k+2 

2.   B  =  (k+l)  B  L 

q    v  -         max   |e  I 
0<i<k-l '  i1 

By  Lemma  2 . h 

|e  |  <  e1"*   n    k_1  (E  +  -iL(t  -  t.  _))  (2.71) 

1  n '  —  1-q   n    k-1' ' 


We  conclude  that  as  h    ->  0,   e   ->  0. 

max      '  n ' 


Q.E.D. 


1*2 

Similarly,   stability  and  convergence   of  the  variable  k-step 
Adams-Bashforth  method  can  be   carried  out   the  same  way.      Thus,   Theorem  2.2 
and  Theorem  2.3  also  hold  for  Adams-Bashforth  method. 

Example  2.3 

Numerical  testings  of  variable   3-step  Adams  method  with   step 
size  being  alternated. 

Again,  machine   computation  is  performed  on  the  same  initial 
value  problem  as   in  Example  2.2.      The   results   are  tabulated  below  in 
Table  2.5   and  Table   2.6.      (One  program  (for  Table  2.5)    is   selected  to  be 
attached  in  Appendix  B). 


TIME 

COMPUTED 

EXACT 

SOLUTION 
69C734~3173~0 

rc~ 

SOLUTION 
~~Z~.  6 90 7 :3"4T3C 60" 

0.370000 

00 

0. 

or 

0.6450CD 

00 

0. 

52466252550 

C3 
00 

0. 52466254210 
0.  39851904H0 

OC 

0.92000D 

00 

0. 

39851902360 

or 

0.1195C0 

01 

c. 

,30270393720 

00 

0. 3C27039542D 

or 

0.  1470C0 

01 

c. 

,22992546950 

0  0 

0. 22992548520 

oc 

0.  174 50  0 

01 

0. 

,1  7464<t97490 

00 

0.  174t.449890D 

00 

0.2C20CO 

01 

c. 

13265545280 

cc 

0. 13265546510 

or 

0.2295C0 

01 

0. 

,10076138270 

CO 

0.10076139330 
0.  76  53  3  54  542  0- 

Z'n' 

C.257CC0 

01 

0. 

,  76535536450- 

■01 

0 .284500 

01 

0. 

58134289210- 

■CI 

0. 58134266740- 

■01 

0.312000 

01 

0. 

,44157162160- 

■01 

0.441571o8420- 

■01 

0.339500 

01 

0. 

33540549010- 

■01 

^.33540584170- 

-01 

0.3670CD 

0  1 

0. 

,25476465710- 

-01 

C. 25476469950- 

■01 

0.3945CO 

01 

0. 

,  19351212920- 

14698641 700- 

•01 

0.  19351216370- 

■01 

0.42  2000 

01 

J. 

■CI 

0. 146986445T0- 

■01 

Q. 449 SCO 

01 

*      i 

,  11  164678350- 

-CI 

0.  11  1^4  6  80  62  0- 

.  .1 1 

0.477000 

01 

0  • 

9480  3783340- 

-C2 

0. d4 80 3 80 1600- 

■02 

0.5C45C0 

CI 

0. 

.64414588980- 

■:?. 

r .644146C364U- 

■02 

Table  2.5.      Integration  of  y'   =  -y   ,  y(Q)  ■  1  using  3-step  Adams 
Variable  Step  Method  with  Step  Size   Alternated  with    .05   and   .005 


U3 


m-r»m  COMPUTED 

TIME 
X  ■*  SOLUTION 

0.370000    OC  0.  167416lPCCD-f)l 

0.6450PD    00 P...1Z3O7680C6P    Op 

0.9200CD    00  0.71639796000    OP 

0.119500    01  0.7C39755401D    01 

0.1A70P0    01  0. 46694833100    01 

0.174500    01  0,92721772510    01 

0.202000    01  0.16649664160    02 

0,229500  .0  1  _ " Q  ,27  7_4  1 5J5  7_3  5  D__02_ 

D.257000  01  0.  43624704010  02 

0.234500  01  0.6  5513240700  P? 

0.312000  01  0.9476P54336D  02 

0.339500  01  0.13284925230  03 

0.367C0D  01  0.18141126720  03 

0.  394  500. ..C.1 0, _2 4 2.2 07 74  7 2D  0_3_ 

0.422000  01  0.31713°1 1060  03 

0.44  9  5  00  01  0_.4C824  303530  03 

0. 477000    01  C. 51769445P4P    03 

0.504500    01  0.64780557660    03 


EXACT 
SOLUTION 

0.  18741610000-01 

-0.t_L13.Q7ft 80 060  oo 

0.7163  92  96  0  00  00 

0.20  39  2  5540  IP  01 

0.466948R31 00  01 

.0,92  72177  2  5  10  01 

0.16649664  160  02 
0.27  74 1 5 5 2  3 5 0 _ 0 2_ 

C.  4 362 4 704 010  02 

.0.65513240700  0*> 

0.9475P543360  02 

P.  .13284925230  03 

0.  18141126720  03 

-Q__24  2  2  0.17  4  72  0  _0  3 

C.?17l"?9ll06D  03 

0,40 874  303 5  3D  03 

0.  51769445R40  03 

_0,  64 780 55 7660  03 


Table  2.6.   Integration  of  y'  =  kt    ,  y(0)  =  0  using  3-step  Adams 
Variable  Step  Method  with  Step  Size  Alternated  with  .05  and  .005 

From  Example  2.2  and  Example  2.3  it  is  seen  that  for  Adams 
methods  the  variable  step  method  is  extremely  stable.   Even  when  subject 
to  adverse  condition  (step  size  being  alternated  constantly)  its  convergence 
behavior  is  superior  to  the  interpolation  method  with  stabilizing 
condition  (step  size  being  held  constant  every  k  steps).   It  should  be 
pointed  out  that  the  programs  in  Example  2.3  are  written  without  explicitly 
taking  condition  (2.53)  into  account.* 


*  Theorem  2.3  was  proved  by  Piotrowski  [12  ] ,  however  his  proof  assumed 


without  proof  that  the  3   .  were  uniformly  bounded. 

n,i 


kk 

3.   STABILITY  MD  CONVERGENCE  OF  GENERAL  MULTISTEP  METHODS 

That  the  variable  step  method  is  more  stable  than  the  interpolation 
method  in  the  case  of  Adams  method  is  apparent  in  Chapter  2.  We  wish  to 
generalize  the  investigation  to  the  case  of  general  multistep  methods.  Two 
independent  proofs  on  convergence  are  given  for  the  interpolation  method 
and  the  variable  step  method,  respectively.   In  section  3-3  numerical 
testings  are  carried  out  for  the  3-step  stiffly  stable  method  in  an  effort 
to  offer  one  more  quantitative  comparison  of  the  interpolation  method  and 
the  variable  step  method. 

3.1  Interpolation  Method 

In  proving  the  convergence  of  Adams  interpolation  method,  we  make 

use  of  three  essential  conditions:   (2.22),  (2.23),  and  (2.2U).  We  observe 

that  not  all  of  these  conditions  remain  valid  when  we  are  dealing  with  a 

general  interpolation  method.  Fortunately,  (2.2*0  still  holds  in  the 

general  case  and  that  along  with  appropriate  step  size  change  constitutes 

sufficient  condition  for  the  convergence  proof. 

Let  h  be  the  maximum  step  size,  i.e.  h  =  h    .   Let  9(t)  be  a 

max 

differentiable  function  in    [0,b]    such  that 

^0   <  A   <_6(t)    <_  1  for  some  A   >   0  (3.1) 

max      |e'(t)|    <  0'        for  some  0'    >   0 
t    [0,b] 

,  h     =  h9(t   ) 
V.  n  n 

Then  from  (2.1+7) 


>+5 


h  e(t. )  -  h  e(t.  , ) 
5i TeTT-^ 


(3.2) 


e'(s.)(t.  -Vl) 

e(ti-l> 


S.     £     (t.         ,    t. ) 

1  1-1  1 


e'(si>  hi-i 


=  h  e'(s. ) 


Thus,   6.    is   uniformly  bounded  by 


6.  |    <_  h0'      ¥  i 


(3.3) 


We  shall  express  the  step  changing  matrix  C.  in  an  alternate  form. 


Lemma  3.1 


The  step  changing  matrix  C.  for  a  (k+l) -value  general  multi- 


value  method  satisfies 


C.  =  I  +  6.  C. 
l        li 


(3.U) 


where 


C.  = 

i 


o 


o 


Z      (f)  5^_1 
1=1     *    X 


Moreover,  C.  is  uniformly  bounded  by  some  constant  k~  >  0- 

1  Vy 


Proof 


O 


U6 


C     = 
i 


o 


(wt)  O 


o         {1+6 . )» 


=   1  +   6 


o 


o 


L 


£-1  * 


J 


1  +   6.    C 
1     1 


By   (3.5) 


'I  X 


k 
E 

k 


k)  is.r1 


<        E      (f)    (h  G'/"1     *i 
~    £=1     * 


hi 


Consequently  for  some   constant  k~  >    0,  we  have 


Cjl    ±k6     V-i  (3.5) 


Q.E.D. 


The  next  lemma,  follows  from  Lemma  3.1. 


Lemma  3 . 2 

Let  S  and  S  be  matrices  which  correspond  to  (2.l6)  and  (2.19) 

based  on  a  (k+l) -value  general  multivalue  method.   Then 

S   C  _  =  S  +  h   .  S 
m  m-1        m-1  m 

where  S  is  uniformly  bounded  by  some  constant  k  >  0. 
m 


Proof 

Analogous  to  (2.18)  it  is  clear  that 

S  =  S  +  h  .  S  (3.6) 

m       m-1  m 

where  S   is  uniformly  bounded  by  some  constant  k~  >  0.   Substituting 
m  S 


(3.M  and  (3.6)  into  S   C   .  we  have 

m     m-1 


S     C     _    =   (S  +  h     _    S    )(I  +   6     _    C     _)  (3.7) 

m     m-1  m-1     m  m-1     m-1  u    " 

=  S  +  hnS+6,SC  +h  6  SC 

m-1     m         m-1  m-1  m-1   °m-l     m     m-1 

h   e'(s        ) 

=  S  +  h     .    S     +  h     ,     Szi_  s  a 

m-1     m         m-1  h     ,  m-1 

m-1 

+  h      -i    <S      .    S      6     , 
m-1     m-1     m     m-1 

e,(Sm    J 

=  S  +  h         S     +  h     ,  ,    m~1.    s  c 

m-1     m         m-1       e(t        )  m-1 

m-1 

+  h     _    5     .  -S-     C     n 

m-1     m-1     m     m-1 

=  S  +  h     ,    S 
m-1     m 


U8 


where 


■  •(■      ,) 

S     =  S     +     Q/^m  \    S   C  +  6      .    S      C     .  (3.8) 

m         m  6(t      , )  m-1  m-1     m     m-1  ' 

m-1 


Thus, 


S  S  +        fl,.lu  'J        S  C     .        +     (5         I        S     C      _ 

m1  '    —   '  '    m1  '         '    8(t     ,  J '     '  '        m-11  '    m-11     ' '    m     m-11 


< 

1    m1 

+ 

6'(s        ) 
i           m-1 

1  QKJ 

m-1 

< 

Is 

1    m1 

+ 

9,(Sm    T} 

i           m-1 

1   e(t    J 

+    1    »/+        \\     llSH     ll6     nH    +    I6      il     llS    M     MC      ill 
1    8(t        ) '    m-11  '         '    m-11     '  '    m1  '     '  '    m-11  ' 

m-l 


ik§  +  T*Vkc  +  h0'  kskc 


Suppose  we  set 


Then 


k  =  ks  +  Xksk6  +  h0*  kskc  (3'9) 


||S  ||  <  k  ¥  m  (3.10) 

Q.E.D. 

Note  that  all  the  matrices  used  here,  such  as  S,  S  ,  C  ,...,  etc.,  vary 

mm 

with  different  multistep  methods.      Hence,   the   associated  hounds   should 
not  he   regarded  as    fixed  hounds   for  all  methods.      We  need  the  ahove  lemmas 
to  arrive   at   the  convergence   theorem. 


Theorem  3.1 

The  general  (k+l)-value  interpolation  method  is  convergent 
provided  condition  (3.1)  is  satisfied. 


^9 


Proof 

At  time  t     the  error  e^  satisfies    (2.lU) 

%  =   SN  CN-1  %-l  +  ^ 
for  the  matrices   S     and  d^  associated  with  the   given  method.      By 
Lemma  3.2 

SH  °N-1  "  S  +  hN-l   §N 

Thus ,  we  have 

e„  =   (S  +  hT  ,    S,T)   e_  .    +  d_T 
-N  N-l     N     -N-l        ^=N 

N 

N-i  ,             -                         x          N 

=     E      S      J  h.    .    S.   e.    .    +  d. )   +   S     p 

,=1  3-1     J  -j-i      -j               -o 

Let  k     and  e   he   associated  hounds   corresponding  to    (2.2U)   and   (2.30). 
Then 

N 

119,11  <  r   !|sH-J  (Vl§.  aa4*4,)||  +  MsHe0||  (3.ii) 

The  remaining  part  of  the  proof  is  trivial.      We   adopt   a  technique  previously 
used  in  Theorem  2.1  to  get 

113,11    K^llsoM    +t)^SkltH"t0>-|  (3.12) 

Suppose  we  use   a  convergent   one  step  method  to   start  up  the  multivalue 

method.      It   is   clear  that   as  h  +0,      |e.T|  I    ->  0   for  all  N  such  that 

max  ' ' — N ' ' 

tN  £    [0,h]. 

Q.E.D. 
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We  like  to  point  out  that  Theorem  2.1  is  not  a  particular  case 
of  Theorem  3.1,  as  the  two  proofs  are  independent  and  distinct.   In 
Theorem  2.1  it  is  essential  to  establish  uniform  bound  for  n  (S  C.   ) 
whereas  in  Theorem  3.1  condition  (3.1)  plays  an  important  role. 


3.2  Variable  Step  Method 

Consider  a  k-step  general  variable  step  method 


£   (a   .  y  ,  +  h„  ,    3_  ,  f „  , )  =  0 

i=0 


n,i  "  n-i    n-1  n,i  n-i 


(3.13) 


of  r-th  order.   We  require  that  the  polynomials  1,  (t  -  t  v)>. 

(t  -  t   )  be  exact  in  (3.13).  For  simplicity  let  us  denote 
n— k 

s    =t    -t  n9  m=0,l,...,  k-1 
n-m    n-m    n-k 

Then  we  have 


(3.1*0 


1  .  ..     1    1 


S    V4.T     0        h     1 

n-k+1       n-1 


r       r      n       ,     r-1 
s    ...   s  ,  ,  .,   0  rh  n  s 
n        n-k+1       n-1  n 


h         h 
n-1       n-1 


r-1 
rh   s         0 
n-1  n-k+1 


—               — -k 

/s 

an,0 

. 

. 

• 

n,k 

§n,0 

. 

. 

n,k 

_         l 

=  0  (3.15) 


According  to  the  algorithm  for  variable  step  method  stated  in  Chapter  1 
(page  11),  we  convert  (3.15)  into  the  following  (r+l)  x  (r+l)  system. 


n-p1       n-p2      ... 


2  2 

s  s 

n-p1       n-p2 


n-p1       n-p2 
^ 


where 


n-1 


.  .     2h      ,  S 


n-1 


2h      „s 


n-1  n-q        "  n-1  n-q 


r-1  r-1 

rh     ,  s  rh     ,  s 

n-1   n-q  n-1  n-q 


M. 
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n»Pn 


n>Po 


n,q. 


n,q 


r+1-2, 


n,0 


n,l 


n,r 


(3.16) 


°   -  Pl   <   P2    < *  '  * <   P£    -  k 


01qi    <    q2    <...<    q^  <.  k 


and  b      .'s   are   polynomials   of  h     ,»•••»  h     ,     and  the   specified  a      .'s   and 

n,i       *  n-1       n-k         *  n,j 


n,j 


's  for  i  =  0,  1,...,  r.   From  (3.1*0 


k 

5        E   h   . 
n-m   .   , .,   n-i 
i=m+l 


(3.17) 


Then  it  is  readily  seen  that 
1  A '     n-k 


QA(6n-l'    6n-2'---'   6n-k+l} 


where   Q     is   a  polynomial  in   6      -,  »    ^      05'.'5    <$      ,  ,-,•      The  existence  of  unique 
A  n— 1        n— <—  n— K+J. 

solutions   for   (3-l6)    in  the  case  of  fixed  step  method  implies   that   Q     is 
nonvanishing  in   a  neighborhood  of   (0,   0,...,    0).      Hence,    for  some   constants   d 
and  D  there  exists   a  compact  region 

"■  {(6n-i---'  W'^Vi 6n-k+i  ^  *}         (3-l8) 


on  which  M  is  non-singular.   This  assures  the  existence  of  unique  solutions 


A 


for  a    , .  . .,  a     ,3     , . .  .  ,  3 

n,p1'     n,p£    n,q1'      n, 


Vi-l-Ji 


in  (3.16). 
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We  shall  show  that  these  coefficients  are  rational  functions  of 


6n-l'"->  6n-k+l' 


Lemma  3.3 


Let  6   .,...,  6  ...  satisfy  (3.18)  ¥n.   Then  the  coefficients 
n-1       n-k+1 

G  „,...,&  ,  »  3  «»■...  p  ,  can  he  represented  in  the  form 
n,0       n,k   n,0       n,k 

<S   .  =  a.  +  R   (6    ,...,  6   .._)  (3.19) 

n,i    l    a.    n-1       n-k+1 

WW  («rt w  (3-20) 

j 

where  R   and  E„   are  rational  functions  of  6   ..,...,  6  ,  ,,  with  the 
a.      g  .  n-1       n-k+1 

i       J 

numerators  containing  no  zero-th  degree  element. 


Proof 

It  is  clear  that 


Pa.  (6n-l'""  Vk+l) 


&     .  =  — -7T ; r  (3.21) 


where  P   and  Q   are  two  polynomials  which  are  relatively  prime.   The 
i       i 

coefficient  a   .  reduces  to  a.  when  the  step  size  is  fixed  throughout  the 
n,i  i 


computation,  i.e., 


a 
n,i 


Pa>  (0,  0,...,  0) 


i 


6   .  =...=  6 •    .  =  0    ai 

n-1        n-k+1 


07  (o,  o,...,o)^tti    (3'22) 
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Apply  partial  fraction  technique  to  get 

P   -  a.  Q 
a .    i  a . 

=  a.  + 


n,i    i       Qa 
i 

In  view  of  (3.22),  P   -a.  Q   is  a  polynomial  which  vanishes  at  (0,  0,...,0) 

i      i 

Hence  P   -  a.  Q   does  not  contain  any  zero-th  degree  element.   Thus, 
a.     l   a. 
l        l 


(3.19)  is  verified  when  we  set 


P   -  a.  Q 
a.    1   a. 

R   =  — ^ (3.23) 

a.       <ci 
1        a. 


Similarly,  we  can  obtain  (.3.20). 


From  (3.2)  one  can  express  6   .  as 

n-i 


Q.E.D. 


6   .  =  h  0 '  ( x   .)   i=l,...,k-l 
n-i         n-i 


n-1  0(t  .  j 


(3.25) 


0»(t   .) 
n-i 

:n-l 

Thus,  (3.23)  and  (3.2*0  yield 

R   (5   .,...,  <5   V4_n  )  =  h    5   . 
a.   n-1       n-k+1     n-1  n,i 
1 

where  for  some  constant  K,.   >  0,  la   .  I  <  K_   ¥n.   Take 

a.       '  n,i'  —  a. 
1  1 

■mo  ■y  I  I 

K-  =  n     . '    n    {K^    }.      Then     a      .is   uniformly  bounded  by  K„, 
a       l<i<k-l        a.  '    n,i '  a 

la      ,1  <  K.  for  all  n  and  i  (3.26) 

1    n,r  —    a 

Combine  (3.19)  and  (3.25)  to  get 

a   .  =  a.  +  h   .  a   .  (3.27) 

n,i    1    n-1  n,i 

Substitute  (3.27)  into  (3.13)  and  employ  a  previously  used 

technique  in  treating  f   .  to  get  (assuming  a  n  =  -l) 

n-i  n,U 


5h 


D  "  (Xiyn-l  ••'  -  Vn-k  =  hn-l  (3n,0  yn  +-"  +  3n,k  ^n-k5 


+  h    A  (3.28) 

n-1  n 

where  3   .  =  a   .+3   .  f  (n   .  ),  n   .an  intermediate  point        (3.29) 
n,i    n,i    n,i  y  n,i   n,i  r  ?/ 

between  0  and  y 

A  =  £    f(0,  t  )  +  ...+  §    f(0,  t   )  (3.30) 

n    n,0       n         n,k       n-k 

It  is  clear  that  6   .  is  uniformly  bounded  by  some  constant  K>  >  0. 
n,i  B 

Let 

K3=Ka+KBL  (3-3l) 

A  =    (k+1)   Kg  L  (3.32) 

Then  Q      .    and  A  are  uniformly  bounded  by  Kn  and  A,  respectively.   The 
n,i      n  3 

above  findings  are  summarized  in 


Lemma  3.^- 

Let  the  step  size  change  be  such  that  (6   ......  6  ,  _  )  e  U 

n-1  n-k+1 

for  all  n.      Then  a  k-step   general  variable  step  method   (3.13)    can  be 

converted  into    (3.28)  with  3      •    and  A     being  uniformly  bounded. 

n  ,i      n 

We  state  two  useful  lemmas  from  Henrici  [9].   The  first  one 
is  the  same  as  Henrici 's. 


Lemma  3 . 5 


k       k— 1 
Let  the  polynomial  p(c)  =  an  ?  +  ai  E    +...+  a  satisfy  the 


condition  of  stability  (namely  the  modulus  of  no  root  of  p(c)  exceeds  1, 

and  the  roots  of  modulus  1  be  simple)  and  let  the  coefficients  y    (£=0,1,2,...) 

be  defined  by 
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1 =  Y  +  Yn  5  +  Y0  C2  +..  .  (3.33) 


a„  +  a,  E  +..  .+  a   C^ 


•o   -1^+-.-+\^ 


Then 


r=    SUP   |v  I  <  « 

£  =q  ,1 . . .  '  h  ' 


The  second  one  is  subject  to  slight  modification.  We  will  offer  the 
justification  for  such  change  in  Lemma  3.7 • 

Lemma  3 . 6 

Consider  the  nonhomogenous  linear  difference  equation 

a_  z  +  an  z  _+...+  a.  z 

0  m    1  m-1        k  m-k 

=  h    {3    z  +...+  g    z  ,  }  +  h  ,  A   (3.3*0 
m-1   m,0  m       Hm,k  m-k     m-1  m 

Suppose  the  polynomial  p(?)  =  a  ?  +...+  a   satisfies  the  condition  of 

stability,  let  B*  ,  3  ,  and  A  be  non-negative  constants  such  that  (for 

m  =  0,  1, .. . ,  N) 

|B  J  +  |8   J  +...+  |3  J  <  B*  (3.35) 

m,0'     m,l         m,k'  — 


m,0'  - 


(3.36) 
A  J  ■   A  (3.37) 


nr 


0<h      -,    3    |a      l~     <q<l     for  some  q   >   0  (3-38) 

—     m-1  0  — 


Then  every  solution  of  (3.3*0  for  which 


satisfies 


zj  1  Z    u  =  0,  1,...,  k-1  (3.39) 


z  |  :  K*  eL*  (tm  "  t0)  (3.U0) 

m 


where 
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kArz  +  TA(t  -  O 

K*  =  — i£ 0-  (3.H1) 

1-q 

L*  =  g-  (3.U2) 

A  =  |«0I  +  \a±\   +---+  l\l  (3.U3) 


Lemma  3 . 7 


Let  z  satisfy 
m 


m-1 


Then 


z  I  <  K*  +  L*  Z     h.   z.  (3. MO 

m1  —         .  _   l  '  i1 
i=0 


z  I  :  K*  eL*  (tm  "  t0)  (3.U5) 

m1  — 


Proof 


We  shall  show  by  mathematical  induction 
1.  m  =  1;  from  (3 .kk) 

\z±\    <_K*  +  L*  hQ  |zQ| 

<_  K*  +  L*  h  K*  by  ( 3  .  Ul ) 


=  K*  (1  +  L*  h  ) 


L*h 
<  K*  e 


=  K*  e 
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2.  Assume  (3.^5)  for  m<_n.  From  (3.UU) 

n 

n        L*  (t,.  -  tn). 
=  K*  +  L*  L     h.  K*  e     X    U 
i=0  x 

n         L*  (t.  -  t) 
=  K*  (1  +  E   L*  h.  e     X    U  ) 
i=0     x 

n  L*  (t,  -  t.)    n   L*  (t.  -  t  ) 

=  K*  (1  +  S   (1  +  L*  h.  )  e     X         U  -  I     e     X    °  ) 
i-0         X  i=0 

n   L*(t.+1  -  t  )    n   L*  (t.  -  t) 
ji  K*  (1  +  £   e  -  E   e  ) 

i=0  i=0 

=  K*  e 
Thus,  the  induction  is  completed  and  the  claim  made  in  (3.1+5)  is  verified. 

Q.E.D. 


We  proceed  to  show  the  stability  of   (3.13) 


Theorem  3-2 


A  k-step  general  variable  step  method   (3.13)   is   stable  provided 
the  step  size   changes   satisfy  condition    (3.l8). 
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Proof 

Let 

B*  =    (k+l)   K3  (3.1+6) 

3   =  K6  (3.U7) 

Y  =       max       ly   I  (3  kB) 

0<i<k-l    'M.1  Kl-io) 

Furthermore,   let   A   satisfy    (3.32).      Choose  T    according  to  Lemma  3.5. 

Given  some   0   <   q  <   1,    choose  h  according  to   (3.38).      Then   it   follows 

from  Lemma  3.6   that   at  time  t 

L*    (tN   "   V 
|yNl    ±K*   e  N  °  (3.U9) 


where 


kAIY  +    TA(t      -   t    ) 
K*  =  . S 2_ 

l-q. 


l*  = 


A  =    |o0|    +    |a1|    +...+    |akl 
We   conclude   from  (3.1+9)    that  the  method  is   stable. 


Q.E.D. 


Finally  we  come  to  grip  with  the  convergence  of  (3.13)  in  the 
following  theorem. 

Theorem  3.3 

A  k-step  general  variable   step  method   (3.13)    is   convergent 
provided  the  step  size   change  satisfies   condition    (3.l8) 
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Proof 

Consider  the  associated  difference  operator  L     for   (3.28). 

n 

Ln    [y(tn>]   =  y(tn)    -  ax  yft^)    ...-o,k  y(tn_k) 

-  Vl   (Sn,0   «**.»   +---+  Sn>*  f^(tn-k!))   +  Vl   \, 

=  0(hr+1)  (3.50) 

max 

since    (3.28)   is  of  r-th  order.      Define  e     =  y     -  y(t    ).      Combine   (3.28) 

n         n  n 

and   (.3.50)   in  a  similar  manner  as  in  Theorem  2.3,  we  have 

e     =  a.,    e     , 
n         1     n-1 

+  ...+  a.    e     .    +  h     .    (|      _   e     +...+  g  e     ,  )   +  h     _    I  (3.51) 

k     n-k         n-1       n,0     n  n,k     n-k  n-1     n 

where  6   .=3   .  f  (c   •  ),  C   .an  intermediate  point  between  y  and  y(t  ) 
n,i    n,i  y  n,i    n,i  *  'n       n 


I    hn-l  .  *  Wl„C«)f.   )  ,*V-1  i  (  v  Ste!)**,,   yt^'fx     ) 


h  ,  k   k  h  /„..  \  _ 

_n=l  v  (   i     n~'T)r  B      Yl   ;(T     ),   x   .  ,  x   .  e  (t  ,  ,  £  ) 

r!   .\\\h  /   n,i-l  y       n,i-r'    n,i'  n,i     n-k'  n' 
x=l  j=l  n-1 

We  set 

i.=  U+1)K6L,       6=K6L,       E  =  0<^-lieil    ((3-52).  (3.53),  (3.5M) 

r     max   i  r +1 1  «        _  /■  o  c  q  \ 

A=Ch    ,pf,,i  y/+J»   for  some  C  >  0  13.55; 

max  te L0,d J   it ; 

Choose  r  according  to  Lemma  3.5-   Given  some  0  <  q  <  1  choose  h^  according 
to  (3.38).  By  Lemma  3.6  we  have,  at  time  t^, 

TB* 

1  N '  —        1-q 

where  A  =  1  +    |  ou  |    +...+    |  ot    |  .      Again,   assuming  the  use  of  convergent 

single   step  methods  to  obtain  initial  points.      Then  it  is   clear  that   as 

h         -»■  0,    lej    ■*  0. 
max  '    N ' 

Q.E.D. 
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It  is  clear  from  the  preceding  theorems  and  lemmas  that  the 
stability  and  convergence  of  the  variable  step  method  depends  on  that  of 
the  corresponding  fixed  step  method.  To  begin  vith  one  must  choose  a 
stable  and  convergent  fixed  step  method.  Conditions  (3.1)  and  (3.l8) 
will  both  be  fulfilled  if  the  ratio  of  successive  step  size  satisfies 
(1.25),  namely, 

h  r-   h   _ 

n    n"-  =  0(h   ) 


h  .,        max 
n-1 

Intuitively,  this  means  that  interpolation  method  and  variable  step 

method  are  both  stable  under  mild  step  size  change. 

The  generalization  of  Theorem  3.3  to  systems  of  differential 

equations  can  be  readily  accomplished  by  replacing  y  with  a  vector  y  , 

absolute  norm  |*  |  with  vector  norm  | |*  | | .   The  final  form  for  error 

bound  will  accordingly  be  the  same  as  (3.56)  with  the  corresponding 

constants  A,  T, .  . .  . 

3.3  Numerical  Testing  of  Stiff  Interpolation  Method  and  Stiff  Variable 
Step  Method 

In  this  section  we  shall  perform  numerical  tests  on  a  multi- 
step  method  other  than  Adams  method.  We  choose  a' stiffly  stable  method 
because  it  is  a  very  useful  numerical  method  in  dealing  with  stiff 
differential  equations.   Furthermore,  the  associated  linear  difference 
equation  of  stiffly  stable  method  is  quite  different  from  that  of  Adams- 
Moulton  method  although  they  are  both  implicit  methods. 

k 

y  =  I   a.  y   .  +  h  R   f  (3.57) 

n    .  .,   1  n-i      0  n 
i=l 

is   fixed  k-step    (k<_6)    stiffly  stable  method 
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(Definition  of  stiff  stability  can  be  found  in  Gear  [8]).   The  corresponding 
interpolation  methods  and  variable  step  methods  for  k  =  2,  3  are  contained 
in  the  following  examples : 

Example  3 . 1 

The  representations  of  interpolation  method  and  variable  step 
method  based  on  the  two-step  Ada.ms-Bashforth  predictor  and  the  tvo-step 
stiffly  stable   corrector. 

(a)      Interpolation  Method 

The   corresponding  matrices.  A  and  i_  of  representation 
(l.lO)-(l.ll)    are 


111 
0  12 
0     0.1 


I    = 


2/3 

1 
1/3 


(b)   Variable  Step  Method 
A-B  Predictor 


*n,(0)  "  Vl  +  Vl  <»£i  Vl  +  ^2  W 


where 


f    .        h    +  2  h 

(p)  _  n-1 n-2 


n,l 


2  h 


n-2 


g(p)  _    Vi 
pn,2     2  h 


n-2 


Stiff  Corrector 


yn,(m+l) 


*S  yn-l  +  «n,l  *n-2   +  Vl  §n!o  '„.<*) 
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where 


a 


a 


(c) 

n,l 

(C) 
n,2 

((C) 

n,0 


n ,  ( m) 


Vl  ^  \-2' 
hn-2   S  Vl  *  hn-2> 

hn-2    (2  Vl  +  hn-2> 

(Vl  +  V2> 
12   Vl  +  \-2> 


Example   3.2 

Numerical  testing  of  Example   3.1 

Machine   computations   for  the  methods   in  Example   3.1  with 
alternating  step  sizes   are   carried  out  on  the  initial  value  problems 
(2.U0)    and   (2.Ul).      The  results   are   summarized  below  in  Table   3.1, 
Table   3.2,   Table   3-3,    and  Table   3 .k. 


TIME 


COMPUTED 
SOLUTION 


EXACT 
SOLUTION 


0.325 
C.6CC 


CCD    CC 
CCC    CC 


0. 

c. 


722 

546 


4313 
6622 


332C    CC 
?<6C    CC 


C. 

C. 


7225273 

5488116 


536C    00 
361C    CC 


0.875 
0.115 


ccc 
ccc 


CO 
CI 


c. 
c. 


416 
316 


69C4 
4623 


718C    CC 
63<=D    CC 


C. 

0. 


<16£62C 

3166367 


1S7C 
694C 


CC 

CC 


0.142 
C.17C 


5CC 

coc 


CI 

CI 


c. 
c. 


0.197 

0.225 


5CC 
CCC 


CI 

CI 


C. 

0. 


24C 

ii£ 

136 

1C5 


3424 

53_2C 
626c 

2824 


S530 
216C 


CC 
CC 


C. 

c. 


24C5C64 
1826635 


16CC    OC 

667D    CC 


0. 

C. 


13  6  7613 

1C53992 


632D  CC 

<A  1C  CC 

122C  CC 

246C  CC 


C.252 
C.28C 


5CC 
CCC 


CI 

CI 


0. 
0. 


799 
6C7 


5846 
2577 


0.3C7 
C.335 


5CC 
COC 


Ci 

01 


0. 
0. 


461 
35C 


C.362 

C.39C 


5CC 

COC 


CI 
CI 


c. 
c. 


266 
gC2 


1917 

219  6 
C1C3 
C256 


5 ier-c  l 

3_5_6  Cz  C  1_ 
75«0-Cl 
C3_5CjHCL 
5C2C-C1 
9  89C-C  1 


C. 

_o_. 

C. 

_c_. 

C. 

c. 


6CC5631 
6C810C6 
*616<62 
35C6435 


27SD-C1 

2_6Ji_C-C_l_ 
636C-C1 
4  ICC- CI 


2664CCC 
2C2A191 


734C-01 
145C-C1 


C.417 

0.445 


5CC 

ccc 


CI 
CI 


0. 
C. 


153 
116 


431P 
5263 


0.472 
C.5CC 


50C  CI 
CCD  01 


0. 


88* 
672 


9761 
1110 


585C-C1 
233C-C1 
366D-C2 
563C-Q2 


0. 


1537519 
1167656 


166C-C1 
697C-C1 


C. 


6  8  7  C  7  1  3 
6737946 


91CC-C2 
999C-C2 


Table   3.1.      Integration  of  y'   =  -y,  y(0)   =   1  using  3-value  stiff 
Interpolation  Method  with  Alternating  Step   Sizes   of   .05   and   .005 
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TIME 

_C._3.2JC  CD 
C.6CCCCC 
0.875CCC 
0.115CCC 
0.14250C 


COMPUTED 
SOLUTION 

CC_         C.7224661712C  OC 

CC  C".=4Ef2691£7D  CC 

CC  C.4167686183C  CC 

CI  C.31654376C2C  00 

CI  C.24C421C5  -<cp  CC 


EXACT 
SOLUTION 

0.7225273536C  CO 

"C ."S4 Veil 636  1C  CO 

C.  4UE62C1S7C  CC 

0. 31663676940  CC 

C. 2<C5CE4632C  CC 


C.17CCCC 
C.197  5CC 
0.225C0C 
C.2525CC 


01  C.1826C44054C    CC 

01  C13E _5____A  E  * C_C^ 

"Ol  C.1C*"22E£9~16D    CC 

ci  c.8CCCfc<3oee4C-ci 


0.  1826825241C  CC 
C.13E7613122C  00 


C.28CCCC 
0.3C750C 


CI     C.6Ci667E19^C-Cl 
CI     C  ,46  1  5353657C-C 1 


C.  1C539S22460  CC 
Q.6CC5e?1275C-Cl 
C.6CE1CC6263C-C1 
0.46iec£2636C-Cl 


V.335C0C  CI     C.35C5449636C-C1 

_0_»162__C  C_  CI     C. 2662456  241 0  ■_ CI 

C.39CCCC  CI     0.2022186587C-C1 

0.41750C „C1 0.1S358E95C3C-C1 

C.445CCC    CI  C.1U653754CC-C1 

C.4725CC    CI  C.886CC76384C-C2 


G.35C843541CC-C1 
C_^2  664SCS734  C  -  OJ. 
C.  2C24  19  1  14  5C-C1 
0.  1537  51516  6C-C1 
C.1U7E56  6V7C-C1 
C.687C713S1CC-C2 


C.5CCCCC  01 


C.6729397971D-C2 


C. 6737C46S99C-C2 


Table 


3.2.      Integration  of  y'    =  -y,  y(0)   =  1  using  2-step  Variable  Step 


Method  with  Alternating  Step  Sizes  of   .05   and   .005 


TIME 
0.3?5000    00 

0.6  loom    0  9 


COMPUTED 
SOLUTION 

0.  117^7927^0-01 
0.1M^7271?in    on 


EXACT 
SOLUTION 

0.11  IS 664 06 20- 01 

0.  12°600OO0O0    00 


0.  H  i^n  ?-n 
0.1 15000 

:  1 

01 

^  .  590  66  77  "U  2  0 
0. 17  56 « 2  26  7 ?0 

00 
01 

n.  5SM8164 060 
■   0. 17490 C6 2 5 00 

00 

01 

0. 142^30 

o .  i  7  n  ■:  *>  o 

01 

ci 

0.41 35501 59RQ 
0.  q*  6^2  794- SO 

01 
0  1 

0.4 123437* 9 in 

o.  e^nccccoo 

01 
01 

0. 197 son 

0.2^009 

01 
01 

0.152  vu«453n 

n .  /  q  a.  r.  n  ?  1 1  s  ^  0 

0  2 

0.  15214P753  90 
-]m  2  5  6  2  a  Q  0  (■•>  2  5  0 

02 

r'2 

0, 2t?600 
0.7  3  0000 


CI 
01 


0.  3  0'7q  ^9 
0  .  "3  *  K  0  )  0 


lV 


O.^^O")    CI 


0. 406*6  8 3  65  20     0? 
0. . h.  1.5.1  ?_6 64 4  jQ    Q2 

7  9 


■-  .  ^z,c  *,*  4  7' 

0.  17'v):  I  97660 


12 
0_3_ 

0.17?  7" 5 07 770    0  3 
-      ?  ^  i  4  i  n  as  >  ?o     r^ 


0.  4QA4fn594l40     C2 

0. 614 S5 600  00  0  Q2 
"> .  ;'  9  4  "  M  a  4  a  i  a.  o  ~  2 
0 ._!  2  69  4  4  5  06  20  0  Q 
0.1726  7602540    03 


0.417500  01 

3.445  Q.,n  Ql 

0.4  775~0  fl 

0.5QijO;)Q  01 


0  .  "*  0  9  n  3 1  r-  5  f-  a  0  0  3 

0.  ^q?2c'°  ^Q9Q0  Q 3^ 

~.4s*t56H~«5in  9-3 

0.625  1^0771 ln  0* 


0.3  02  9.2  6  6*790  0  9 

Q. J 9~>1  9 900 6 2 0  03 

n. 49R43  351169  0  2 

0. 6250C0000Q0  03 


Table   3.3.      Integration  of  y'   =   Ut    ,  y(0)   =  0  using  3-value  Stiff 


Interpolation  Method  with    Alternating  Step  Sizes  of   .05   and   .005 


6U 


TIME 
0 . 3? ^r ^n    no 

o.6  00ooo~~ooT 

0.  H7^nv)n    00 
01 


COMPUTED 
SOLUTION 

).Po^7,ri^in  oo 

_n.j_H'*^7i^^^  oo 

0.  L  7  «^ 61^7^0  0! 

0.4  I.'  3S?1  15  70  0] 


EXACT 
SOLUTION 

_0,  1  1  lS6ft4C6?0-01 
0.  1 2 9600000 O^  00 
0.  sqm  P]  44050    np 

0.  17^''0C^?c-0^  "ci 

0,  <*  1?3^^7H^10     01 


0.17  0''  10  01 

0.197*00  01 

o.??5co*n  01 

Oj25?^lD  01 


o.  ^?f.n7^Pi^n    oi 

0.  lr^?6S345?0     0? 


0,  H  1^21  VK  *n    ."1 
0.  15?14fl75?9n    0? 


O.?«0000    01 
0.3 07* on    01 


0.25*441  175 
0.4C667763? 
0.61 4*9  1  $44 

0.  fH}4"*7301  4 


^n  o?  o.?e?6?»")0625n    o? 

RO     02  0.  A064R5941 '.0    0? 


AO     0?  O.M44550000O     0? 

in    Q2  0. 39403344140    02 


0  .  3  3  5  ~  ^  n 

n 

0  .  1  2  ^  7P> ")°  ^Qn 

03 

0.  12  5Q445C620 

0  3 

0. 3  6?K on 

01 

0  .  1  7  ?  7 1  5  6  0  0  n  n 

03 

0. 17267602*40 

0* 

o.**Qoooo 

01 

0.  ?3  1  ">P':o  ?ocn 

°3 

r,21l  i44100n(i 

03 

0.4175 o  n 

01 

0.303R7'321  i  an 

0^ 

0.  303.a26A.R79n 

0^ 

0. 445C0n 

01 

0.  *c2l9R6O06n 

03 

0.  332139Q0'-?n 

0^ 

0 .  4  7  ^  5  r  n 

r\ 

0.4QPS^OP^4rn 

03 

o. 49«43^531 60 

r3 

0.50000n    01  0. 626075^7 7°n    03 

Table   3.k.      Integration  of  y '    =   kt    ,  y(0)   = 


0. 6?50C00C03O    03 
0  using  2-step  Stiff  Variable 


Step  Method  with  Alternating  Step  Sizes   of   .05  and    .005 

Example  3.3 

The  representations  of  interpolation  method  and  variable  step 
method  based  on  the  three-step  Adams-Bashforth  predictor  and  the  three- 
step  stiffly  stable   corrector. 

(a)      Interpolation  Method 

The  corresponding  matrices  A  and  l_  of  representation 
(l.lO)-(l.ll)    are 


A  = 


1111 
0  12  3 
0  0  13 
0     0     0     1 


L  = 


6/11 

1 

6/ll 
l/ii 
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(b)   Variable  Step  Method 
A-B  Predictor 

Already  presented  in  Example  1.2 

Stiff  Corrector 

For  simplicity  we  shall  determine  the  coefficients  in  terms 
of  time  intervals. 

y    /  +1x  =  a(Ci  y   i  +  «(Co  y    o^Ccl  y    »  +  h   i  &{cl  f    t  \ 

*'n,(m+l)    n,l  n-1    n,2  n-2  n,3  n-3    n-1  n,0  n,(m) 


where 


a(Cn}  =  (l/|M.|)(t  -t    )(t  -t   J2(t   -t    )[(t  0-t      J(2t  -t   _-t   _) 
n,l      'A'   n  n-1   n  n-3    n-2  n-3    n-2  n-3    n  n-2  n-3 


"  <W3)2] 


Ac) 


a  '   =  (l/lM.  )(t  -t   .  )(t  -t   J  (t  _-t   ,)[-(t  ..-t   _)(2t  -t   -t   _) 
n,2      '  A1   n  n-1   n  n-3    n-2  n-3     n-1  n-3    n  n-1  n-3 

n  n-3 

-,    .(C)    .(C) 
n,3       n,l    n,2 

HI  -   (1/lMAl>(Vl-V3)(V2-V3)(VV3)[(V2-V3)(tn-V3)(VV2) 

"  (tn-l-tn-3)(tn-tn-3)(Vtn-l)  +  (tn-l-tn-3)(tn-2-tn-3)  ^n^n-l5  ] 

|M.  |  =  (t  -t  . )(t  _-t  ,)(t   -t   _)[(t  0-t   _)(t  -t   ,)(3t  -2t   -t   ,) 
'A1     n  n-1   n-1  n-3   n-2  n-3    n-2  n-3   n  n-3    n   n-2  n-3 

-  (t   -t  _)(t  -t  _)(3t  -2t   -t  J  +  (t  _-t  _)(t  0-t  _) 
n-1  n-3   n  n-3    n   n-1  n-3      n-1  n-3   n-2  n-3 


(-Vl+tn-2>] 


f       =  f(y   .  .,  t  ) 

n,(m)     'n,(m)   n 
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Example   3.^ 

Numerical  testing  of  Example   3.3 

Again,  machine  computations   for  the  methods   in  Example   3.3  with 
step  size  "being  constantly  alternated  are   carried  out  on   (2.U0)   and  (2.Ul). 
The  results  are   summarized  in  Table  3.5,   Table   3.6,   Table   3.7,   and  Table   3.8, 


COMPUTED  EXACT 

SOLUTION  SOLUTION 


TIME 
0.325COO    00  0.72252631980  0. 72252735  ISC    00 


0.600C00    0^  C.  54F8556493H    00  C. 54 88 1  16 36  10  00 

(K875C00    00  0.41546214340    00 0. 4 1 68 62 C  1  9 70  CC 

0.115C00    01  6. 3616<52<?0"C3D    00  0.31663676940  CO" 

0.  1  4  2  5  0  0_  01 -  C  .  12C9  0_5  C4840    CI  0,  2  4C5C846370  CC 

0.17 CO 00    01  0. 46 813  9 29  160    0  2  0.  1826  83  52410  CC 

0.1975CO    CI         -0.  150C ?7C?56D    C4  0.  1  3876  1  3  122C  CO 


0.225000  01  0.48272478480    05  C. 1C 53 9 9 22460    CO 

0.252500  01  -0. 15530525940    07 0  .  P005 8 3  1  2 790-C 1_ 

0~.28CC0O"01  C.49~96590096"6    C8~         0. 6C810062630-6i 

0.107500  01         _-0.16CJ_5^QC4CO_  1C 0.46  18962  F3  PO-C  1 

"6. 335C0C !  01~  0V5 17 188  42  330    11  0.35084354100-01" 

_?_»_Lft  25rO  01  -0.  16639*49070     13  0.  2  6  64  qp  C7340-0  1 

0.390000  01  0.53533781770    14  0.20241911450-01 

0.4175C0  CI -0.17223103220    16 0_._15  375  1  9  166C-01    _ 

0.445C00  01  0.5.541.1376750    17  0.  1 1  678  566970-C1 

_0._4725C0  CI -0, 17_8?73?7_120_19 0.,  F  F  707  1  39  10D-C2 

C .  5  C 10 C 00  0  1         '"  C.'  57 35 5  3  6 C 47 o"    2  C  3 .  67 3 7 946 99 90 -C 2 

Table  3.5.      Integration  of  y'   =  -y,  y(0)   =  1  using  U-value  Stiff 
Interpolation  Method  with  Alternating  Step  Sizes   of    .05   and   .005 


TIME 


c.37cccr  cc 
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COMPUTED 
SOLUTION 

c.e<?c7?c;667er  cc 


EXACT 
SOLUTION 

0.6<5C73*??C6C    CC 


0.6*5CCD  CC     C.  c2*66*?PCcO  OC     C.ei2*662S*2 

o . 9 2oncr  co         c^cp57Cciccp   cc         c.  2<;f51<;c*i 


0.1  10 
0.  1*7 


sec   CI 
cc^ci 


1C 
in 


CO 

cc 


c. 

r  . 


?C27C57Slin    OC 

7  ?  c  c  ?  "i  1  c  7  2  C     r  C 


C. 

r  . 


2C27C 
22<^2 


355* 
5*85 


2C 

20 


CC 

CC 


0.17*50C    01  C.  17*6*652<5c;0    CC  C. 

0.2C2Cpn    CI  .C.  l?7fe56PlP3r    OC  c. 


1  "5464 
1  3  2 6  5 


CD    CO 

IC   cc 


0.229 
0.257 


ccc   r i 

CCC    CI 


0. 


\CClt2cc.ccX    CC 
76  536537CCn-Cl 


C. 


1CC76 
76c35 


I  3  c  3 
5*?* 


3C    CC 
20-C1 


C.2P* 

0.31? 


5CC    CI 

rrr,    f\ 


C. 

r  . 


5613^5  555 
**1  C7P62  16 


r-ci 

n-ri 


0. 


5  p  1  2  * 
*.  *  1  5  7 


2667 
16P* 


*C-C1 
2C-C1 


0,3?9cor    CI  C.335*l!?732r-01  0. 

C.'i67CCC    CI  0.2**76^^3^-01  Q. 


3  3  5  AC 

2  5*76 


55*1 
*6<=5 


70-C1 
5C-C1 


0.394  500    CI 

c.*?2ccr   c 1 


C.  lc3c16CC36n-Cl  C. 

o.i*ftq»ct^6*2r-oi         c. 


is^5i 
i*6<;e 


2163 
6**c 


7C-C1 
CC-C1 


C.**9 
Q.*77 


5CC 
COJL 


CI  C.  111649^2920-Cl 

01  C.fi*fC5e'a*P?n-C2 


C. 


1116* 
F*K2 


6FC6 
6  C  1  6 


2  C  -  C  1 
CC-C2 


C.5C*5CD    CI 


0.6**1623675f-C2 


0.6**1*6C36*C-C2 


Table  3.6.      Integration  of  y1    =  -y,  y(0)   =  1  using  3-step  Stiff 
Variable   Step  Method  with  Alternating  Step  Sizes   of    .05   and    .005 


TIME 


C.32 
C.6.C 


5C0O 

cc  cr 


cc 

rr 


0.  11 
".11 


COMPUTED 

EXACT 

SOLUTION 

SOLUTION 

134520920-01 

0.111566*C62C-C1 

^ccc7A?;n    "c 

0.  1296.rcH00r>p    CC 

0.R75CCC    CC  0.5597P52216C  CO  C.  5  5 6  1  f  16*C6C  CC 

0.115CCC    CI  C.  2*^*22*1 CCC  CI  0 . 1 7*^0 C62 5CC  CI 

0.1425 CO    rl  -C.!6F=M'75cn  C2  C*123*°7f(nn  CI 

C.17CCCC    CI  C.60C70'a6-g31C  C  *  O.P3c2  1CCCCCC  CI 


C.19 
0.2? 


7CCD 

cCQr 


CI  -0.1666571 5160    r  5  C. 

ci         c  .*7CPg;6*c7^n   C6         c. 


15 


21*F 
62£5 


753S 
C  62  5 


C    C? 
n    C2 


0.25 

0.?* 


25C0 

rr  oo 


CI 
CI 


•0.  13 
f.37 


27c8?^6*^    CP  C. *0**PcS*  1*C    C2 

*3~'a7  6''2n    rc:  C.614656CCC0C    C2 


0.^0 
C.33 


7^or 

5  CCO 

fi.3625^0 

c.3qcccr 


CI 

ri 

X_L 


- 0  •  10557226 P  3  C 
C.  2C75*C 73*PC 

.  r  #   p  "3  p  p  c  /j  C  A  T  «=  0 


11 
12 


C. 
0. 


PS 
12 


I  3 
15 


0. 


17 


*CPP 
55** 
2  6  7  6 
13** 


**  1* 
5C62 

p25* 
1CCC 


0    C2 

C_C3_ 

r    £3 

C     C? 


C.*l 
0.** 
0.*7 

C.5C 


7^cr 

5mn 


CI 

C  L 


0.6  6^C5267CC0     1ft  0. 

0.1PPCl*'a7?2n     IP  c . 


■a  r 

—   ^ 

-  c 


2526 
7135 


6P7C 
CC62 


r    c? 
C    C3 


7CCC 
CCCO 


CI 
CI 


CI* 


CC9'»1635r     IS  0.*SF*225?16 

§  *  5600  KH    21  C.62^CCCr0C" 


0    c? 


Table  3.7.      Integration  of  y'   =  ^t3 ,  y(0)    =  0  using  U-value  Stiff 
Interpolation  Method  with  Alternating  Step  Sizes   of    .05   and    .005 


TIME 

C.37CCCP 
5  CC  C 

cccc 


rr 


C.6A 
C.<3? 


CC 
CC 


0.11^00 
0. 1A7C0P 


CI 

01 
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COMPUTED 
SOLUTION 

0.167F  7U<?f  5P-C1 


0.17 
0.  71 


3  1 5  4  9 

f  CC*6 


Cl^O 

C3in 


00 

cr 


EXACT 
SOLUTION 

C.  1F74U1CCCC-C1 

^ecctr  cc 
2<;*rcr  cc 


c. 


17?C1 

7U3C 


0.  *c 

0.^6 


■3C  -jq  f. 


5F60 
•53  7H 


CI 
CI 


C. 

0. 


2  C  ?  c.  2 
466-C4 


C5A  ( 

FFF  1 


IP  CI 
CC  CI 


c .  1 7  a  c  r  c 

0.202C00 


rl 
01 


C.92 
0.16 


7 ;  ?  e c 

Mccc 


4<57n 


CI 
C2 


C. 

c. 


c2721 


772? 
66M 


1C 
60 


r  i 
C2 


C.22 

0.25 


9C,CC 

7ffcn 


CI 
C! 


0.2  7 
C.  4? 


741  *2 
625CC 


S7C0 

^ccn 


C2 
02 


0. 


277415521 

436 247C4T 


5C  C2 

in  c? 


o.?"^cn 

C.312CCn 

CI 

CI 

0  .  t  5  c  1  a  «5  7  c  i  u  n 
C.9475GC!4^0 

C2 

C2 

C. 6551?24C7CO 
C  .  c  4  7  5  F  5  A  3  ?  6  C 

C2 

C2 

C.  33^50r 
0.367CCT 

CI 
CI 

^.1  ?2F4q655F0 
0.18!41170^?r. 

C^ 

C.  132P4S?C2?C 

c  iei4ii2672r 

C^ 
C" 

C.39A50T 
0.A32C0P 

CI 

01 

C.2422CF215PC 
0. 2171'9f 117Q 

C3 
C3 

C. 2422C77472C 
C.2l71?cnc6D 

C3 
C3 

0.44 
C.47 


9*>C0 
7CC0 


Cl 
CI 


0.40 
C,  51 


F  2  '♦  "*  5 

76C«=^ 


AR9C 


03 
C3 


0. 

r 


4CF24 
517*c 


3C35 

445* 


3C 
AC 


C3 

r  ? 


0.50ASOO  Cl 


0.647FC617540  C3 


C.647fiC5576fn  C2 


Table  3.8.   Integration  of  y'  =  Ht  ,  y(0)  =  0  using  3-step  Stiff 
Variable  Step  Method  with  Alternating  Step  Sizes  of  .05  and  .005 


Programs  for  Table  3.5  and  Table  3.6  are  attached  in  Appendix  C 
and  Appendix  D,  respectively.   Based  on  the  above  numerical  tests,  the 
stiff  variable  step  method  "seems  to  be"  more  stable  than  the  stiff 
interpolation  method.  However,  one  should  be  cautioned  that  this  is  by 
no  means  a  quantitative  conclusion.   This  writer  feels  that  further 
theoretical  research  along  this  direction  may  eventually  lead  to  a  clear- 
cut  answer,  one  way  or  the  other,  such  as  in  the  case  of  Adams  method,  to 
the  relative  merits  of  stiff  interpolation  methods  versus  stiff  variable 
step  methods . 
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k.      PRACTICAL  CONSIDERATIONS  OF  THE  STABILITY  AND  CONVERGENCE  THEOREMS 

Step  size  varying  is  not  merely  of  theoretical  interest  as  further 
understanding  of  the  mult i step  methods  is  enhanced  hy  a  knowledge  of  the 
stability  criteria.   Even  more  significant,  it  is  of  practical  use  in  that 
computer  time  is  saved  hy  reducing  the  number  of  steps  "while  still  achieving 
a  desirable  accuracy. 

Consider  the  error  bound  (3.12).  We  write 

I  KM  i4x)  +  Ki2Z)  e  {h-1] 

where 

4I}  =  KS  Me0H  {h-2) 

t  t  \         -i  K  K    ( t.T   —  t  -  ) 

K(D=I(l+eB  I  0,  {k3) 

*-  K 

Although  the  constants  K    and  K    are  clearly  independent  of  the  step 
sizes  taken,  they  are  not  computable  for  most  problems.  Nevertheless,  e 
is  closely  associated  with  local  truncation  error.   In  practice,  we  shall 
make  use  of  this  latter  property  to  control  the  global  error  cl  . 

Assuming  the  exact  solution  is  sufficiently  dif ferentiable, 
there  exists  a  smooth  function  <j>_(t)  called  the  principle  error  function 
such  that  the  local  truncation  error  satisfies 

^  -   ±<V  CI   +  °(hn-l>  (U-U> 

for  a  r-th  order  method.   Given  some  specified  e  we  choose  the  largest 


possible  h  .  such  that 
n-1 


|d  |  z   h    e 
1  — n '  '  —  n-1 


according  to  (2.30).  We  shall  show  that  such  a  choice  for  h  ,  does  not 

n-1 

violate  condition  (1.25). 
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From  ( U . U )   we  have 

lli(tntl)||  hnr+1  +  o(^+2)1£hn  (U.6) 

If  we   choose  the  largest  h      and  h  so   that   these   are  equalities, 

h     -h     .        .      ||l(t  +.)||   hr+1  +  0(hr+2)   -   H^t   )||    h^1  +  0(hr+2) 
n  n-1    .  1      '  '  *-    n+1    '  '      n n MI-    n    '  '      n-1 n-1 

h         " "" e  h     . 

n-1  n-1 

•  -(I l±{t   *i)||    »     ^  +   0(«     hr+1)    -    I  liK*    )|l    ^   i    +   0(hr+^)) 

E      "-1-    n+1    Mnn  nn  '  '  *-    n   '  '      n-1  n-1 

:    0(h        ) 
—         max 

Therefore,  when  the  local  error  is   controlled,   the  step  size 
ratios    satisfy  the   condition   for  stability  and  Theorem  3.1   shows   that   the 
global  error  is  bounded  and  converges   as  e    goes   to   zero   according  to 
(U.l)    and   (k.2) 

The   same   thing  can  be   said  for  the  variable  method  since  the 
error  bound   (3.56)    can  be   converted  into  the  same   form  as    (U.l).      The 
step  size   changing  mechanism  is  also   similar,  which  we  will  not  repeat  here 
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APPENDIX  A 
PROGRAM  FOR   THE   INTEGRATION   OF  y'    =   -y ,   y(O)    =   1 
USING  U-VALUE  ADAMS   INTERPOLATION   METHOD 
WITH  STEP  SIZE  BEING  ALTERNATED  ONCE  EVERY   THREE   STEPS 


73 
A-VALUE  ADAMS  INTERPOLATION  METHCn  FCR  CY/CT=-Y 


C   STEP  SIZE  INVERTEC  EVERY  THREE  STEPS 
IMPLICIT  REAL*B(  A-H,C-Z  ) 


DIMENSICN  AP<A) , AC(A>  ,ACC(A) ,6<A) ,C(A) 
_f_CY,T_i5-Y_ 

KCLNT=0 
H=0.05DC 


IS  THE  TOTAL  NUMBER  OF  STEPS 
N  =  2 


T=2.0C*H 
ALPHA=0.1CO 


ESTABLISH  THE  NCPSIECK  VECTOR  AT  THE  INITIAL  POINTS 

ACC( 1 )=CEXP(-0.1CO) 

ACC(2)=F*(-CEXP(-C.1C0  )  ) 

ACC( 3)=C.5*H*H*DEXP<-C. 1001 

ACC ( A )=  I [  1 /6  . "DO ) *H *H*H* ( -CEXP ( - C • ICO ) ) 
PRINT  900 


SOC   F0RMAT(//7X, 1HT  15X,2HYN  18X,7HEXACT  Y  15X) 
C   B  IS  THE  NCRCSIECK  VECTOR  AFTER  STEP  CHANGE 
100   KGLNT=KCLNT+1 
ALP=1.00 


IF(KCUNT.EC.l)  ALF=ALPHA 
B( 1)=ACC<1) 


B(2)=ALP*ACC(2) 
P(3)=ALP*ALP*ACC(3) 


R(A)=ALP*ALF*ALP*ACC(A) 
AP  IS  THE  NORDIECK  VECTOR  AFTER  PREDICTION 


ap(d  =  b(i  )  +  e(2  )  +  em*e(4) 

AP(2)=  8(2)+2.DC*B(3)+3.DC*B(A) 


AP(3)  =  B(3)*3.C0*E(  A) 
AP(A)=e(A) 


H=H*ALP 
T  =  T+H 


C(I)  IS  THE  AMCUNT  TO  BE  CCRRECTEC  IN  THE  CORRECTOR 
G  =  H-»F(  AP(  l_)_i  TJ  -  A  P  (  2  ) 


C(1)=(3.C0/8.C0)*G 
C(2)=G 


C(2)=( 3.CC/A.CC)*G 
C(A)=(1 .00/6.C0)*G 


AC  IS  THE  NGROSIECK  VECTCP  AFTER  FIRST  CCRPECTICN 
__  CO  11C  1=1,  A 
11C   AC (I  )  =  £P(  I  )+C(  I) 
GG=H*F( AC< 1) tT)-AC(2) 


C(l)=(3.CO/8.CO)*GG 

_C(2)=GG 

C( ?)=(3.DC/A.CC)*GG 
C(A)=(1.C0/6.C0)*GG 


C   ACC  IS  THE  NCRCSIECK  VECTOR  AFTER  SECCNC  CCPRECTICN 

CO  120  J=1,A 

120   ACC(  J)  =  AC(J  MC(J  ) 

C__JCEANGE      THE    STEP    SIZE    CNCE    EVERY    3_  STEPS 

IF(KCLNT.LT.3)  GO  TO  1A0 
ALPFA=1  .CO/ALPHA 


Ik 


KCUNT=Q 


C   RECORD  THE  INTEGRATION  HISTORY  EVERY  10  STEPS  ALONG  THE  WAY 

_1A0   IF((N/iO,00-N/10).GT.O,0100)  GC_JQ_L5_0 

YN  =  ACC(  II 

YEXACT  =  OE>P<-T) 


PRINT  920,TfYN, YEXACT 
Q?n   FrPPAT<D12,5.2C2C.10l 


15C   IF( T.GT.5.DC)  GO  TC  200 
N_*  N  +1 

GO  TC  ICO 
_2Q.Q._£0NI_IUUE 

STCP 

END 
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APPENDIX  B 

PROGRAM  FOR  THE  INTEGRATION   OF  y'    =    -y,   y(O)    =   1 

USING  3-STEP  ADAMS   VARIABLE  STEP  METHOD 

WITH  STEP  SIZE  BEING  ALTERNATED   CONSTANTLY 
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C   THQrp  STfP  VAR  IARLF-STEP  A-R-M  METHOD  FOP  FIRST  OROFP  EOUATION 


_F  ( 
H  = 

Y! 

Y? 


PL1C  I 
Y*T)  = 

0.05D 
=  1  .no 
=ne  xp 
=  DEXP 


T  RF  AL -8(  A-H,0-Z  ) 
-Y_ 

0 


(-0. 
(-0. 


05001 

1000) 


c  m 

_C T 

C   A 


IS 

LPHA 

N  = 
T  = 
AL 


THF  T 
TM^EjC 
IS  T 
2 

?.  no* 

pha  =  o 


OTAL  NUMRF.R  OF  STFPS 

I1RRENT  T  I  ^F 

HF  STEP  C H  A N G  I  N G i " F  AC TOR- 


H 

.inr 


FT  H 
NTFP 


MltHM2»HN?  HF  THRFF  CONSECUTIVE  STEPS  OVER  THF  TIME 
VAl  N-lfN-2fN-3 


FT  R 
HN 
HN 

MM 


PI,  RP 

1  =  0.0 
2=0,0 

i=o.  n 


2,RPR  RE  THF  COEFFICIENT  CF  A-B  PREDICTCR 

5  0  0_ 

500 
50  0 


- 


FT  Y 
NTFR 


130 


YM 
Y_N 
YM 
pp 
RP 
RP 


N  ,  Y  N  1 
VAL  Nl 

1  =  Y? 

2  =  Y1 


,YN2,YN3  RE  THE 
,N-1 fN-2» N-3 


VALUE  OF  Y  OVER  THE  TI^E 


R,=  YO 

1  =  1.00+(HM1«(  ?.  DO  <HN1+6.D0*:HN2  +  3.D0*HN3  )  )  /  (  6  .  D0*HN2*  (  HN2  +  HN3  I 


C   C 


T  = 

OMPtJ 


2=-  (H 
3  =  (HN 

T+HN1 
TE  TH 


Nl~(  2.00^HNl+-3.D0i-HN2+3.  00-HM3  )  »  /(  6  .  00  *HN2*HN3  ) 
1M  2.r.0'-HNl+3.C0*HN2  I  )  /  (  6  .  D0*HN3*  (  HN2  +  HN3  )  ) 


C   L 
C   I 


ET  Y 
NTFP 


N,YM1 
VALS 


F  PRECICTED  VALUE  OF  Y 

, Y l\> 2 , . Y N 3 "  R E  f HE  VALUES  OF  Y 

N,N-1  ,  N-2,~~N-3 


OVER    THE    TIME 


YN=YN1*mni  ■*<  hpi  -F  (  YN1, T-HNl  )+RP2*F(YN2,  T-HN1-HN2) 
+  «P3,cf  (YN3,  T-HNl -HN2-HN3  )  > 


C       LET     PCS, PC1,RC?,RC3    RF    THE    COEF F I  EC  I FNTS    CF-  A-M    CORRECTOR 
A=HM2  /HNll 


R  = 
RC 


( HN2+HN3 ) 
3=( l.PO+2 


/HN1 
.00--  A 


)/(12.00~B-(i.D0+B)*<B-A) ) 


RC 
RC 


2=-(2.005 

1=0. 50  0- R 


0=1 .no-RC 
r t    the   v  a 


RC 
0  P  D  r- 
~YN=YN  1  +  MM  J. 
C 


R  *■  1  .  0 
C  3  -  (  1 
^-Rf  ? 
I  IJF.  0 
(3C0-* 
+  RC  ? 


0)/(12.00vAMl.D0  +  A)*(B-A)) 
•  ■? O-tiL'  " BC  2--M  1.D0+A1 

-RC1  "'"" 

F    Y    T  k  I  C  E 


F( YN,  T  I+RC1»F ( YN1 , T-HNl J 

~F  (YN2  ,  T-HN  1-HM2)  4-Rr.3*F(  YN  3  ,  T-HN  1  -H  N2-HN  3  )  ) 


C       R 


YM=VM1+HN1* 

r 

EfJnR "~ 
_  If 
YE 

P 


r.    THF     INT 
(  (  M/IQ.PO 

:xact=oexp 

'PINT    020, 


F ( YN, T)+ PC 1«F (VN1 , T-HNl) 

'•F  (YN2,T-HN1-HM2)  »BC3*F  (  YN3  ,  T-HN1-HN2-HN3  )  1 

ITM  history"  every  10  steps  along" the  way 

I.GT.O.OinO)  GC  TO  150 


{  RCO'. 
_+RC2 

egpat 

-N/10 
(-T) 

T,Y*',YEXACT 


150 


FQ 

IF 


RMAT( 012. 
( T.GT.5.0 


N  IT  I 

HN 


ALI7E  THF 
*=HN2 


5,2020. 10) 
0_)  _G0  _TD  200 
DATA  NECFS'SARY 


FOR  THE  NEXT  STEP 
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HN?=HN1 

YN3-YM? 

HN1 

YM2=YN1 
YN1=YW 

C         CHANGE     THE     ST 

C  D 

S  I  7  c 

ALPHA=1 .no/ ALPHA 

GO  tr  no 

?00       CONTINUE 
STHP 

FNO 

APPENDIX  C 

PROGRAM  FOR  THE  INTEGRATION   OF  y'    =   -y,   y(0)    =    1 

USING  U-VALUE   STIFF  INTERPOLATION  METHOD 

WITH  STEP   SIZE  BEING  ALTERNATED  CONSTANTLY 
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C       4-VALUE    STIFF     I  NT F P FOl  AT  ION    METhCC    FOR    FIRST    ORCER    EQUATION 
IMPLICIT    R£AL««( A-HtO-7)  

"  D  I>  E  N  S  I  r  N  'A  P  (  4")  ,  S  C  (4 )  ,  S  C  C  (  4  ) ,  R  (A  j  ,"  C  <  4) 
F(Y,T)=-V  __ 

~  H=C.0  5CO 
C   N  TS  THE  TOTAL  NUMBER  OF  STEPS 


N  =  2 
T=2.0C*H 


ALPHA  IS  THE  STEP  CHANGING  FACTOR 

AlPHA=C.in^  

ESfAPLTSF  t>F    NdRSrECK    VECTOR    AT"  THE     INITIAL    PCTFT 

SCC( 1 )=PcXF(-0.1CQ) 

scc(  ?)-h^(-dexp(  -c.ino  ) 

scc(i)=o.5^F^H^nExp(-c.ino) 


SCC(4)=(  1/6. OC  |*H*H*H*  (-CEXP(-C.IDO)  ) 

_PRINT  °00  

~9"00   FC"R v  ATT/77 X  ,"1  FT  15 X  ,  2>YN  \WX  ,  7"FEXTCT'~Y  15X> 
C   R  IS  THE  NHRPSIECK  VECTfR  AFTER  STFP  CHANGE 
100   P(1)=SCC(1) 

B(2)=ALPHA*SCC(2) 


B(?)=ALPHA*ALPHA«-SCC(3) 
P(4)sALPI-A*ALPHA*ALPHA+SCC(4) 


AP  IS  THE  NORCIECK  VECTOR  AFTER  PPECICTION 

ap(  i)  =r(  i  )*p(  2i*ei  ^i*ehi 


AP(2)=  P(?)*2.CO*P(3W3.ro«P(4) 
AP(3)=B ( 3 )+3.rC-P(4) 
AP(4)=B(4T 

H=H*ALPHA 


T=T  +  H 
C(I)  IS  TFE  AMOUNT  TO  RE  CCRRECTEC  IN  THE  CCRRECTCR 


G=H*F( AP( 1  )  ,T)-AP<2  ) 

C(  l  )  =  (6.nc/n _ c r-)*c 

C(2)=G 

C(3)=(*. 00/11. CO  )*G 


C(4)=(  1  .DO/11. DO)*G 
SC  IS  TFF  NCRCSIECK  VECTOR  AFTER  FIRST  CORRECTION 
HO  110  1=1,4 

11  Q_  SC(I  )  =  AP(  I  )  +C  (  I  )        ___ 

GC  =  H*F ( SCTT)  ,  T  f-SC  (  2  ) 

_CJ1)  =  (6.DC/11.C0)*CG 

C(  2  )  =  GG  "" 

C(3)*U  .cc/u  .no  >*co 


C(4)  =  ( i .nc/n.oo  )  *cg 

C   SCC  I_SJFE  NHR.DSIECK  VECTOR  AFTER  SECONC  CCRPECTICN 
DC  12  0  J  =  l  ,4 
12C   SCC(_-SC(  J  UC(  J)  

r     change"    t hf  "step  size 

ALPHA=1 .OO/ALFHA 


C   PFCOPO  THE  INTEGPATICN  HISTORY  EVFRY  10  STEPS 

IF(  (N/IO.DO-N/IO.OT.O.CICO)  CCJTO  15C 

~YN  =  SCC(iT  ">" 

JT  E  X  A  CT=  C E  X  P  <  -  T  )_  _____ 

PRINT  920ttf YNt'YEX ACT* 
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920 

FDRWAT(C1^. 5, 2020.10) 

150 

IF<T.GT.5.C0)  GC  TC  200 

GC  TT  100 

20C 

CCNTIMJP 

STOP 

FNO 

APPENDIX  D 
PROGRAM  FOR   THE  INTEGRATION   OF  y'    =   -y ,   y(o)    =   1 
USING  3-STEP  STIFF  VARIABLE  STEP  METHOD 
WITH  STEP  SIZE  BEING  ALTERNATED  CONSTANTLY 
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c   JFRCF  S7FP  STTPF  V  A  P  I  AR  L  F-  S  T  F  g  MF7HCD  FOR  FIRST  CROFF  FCLATiriN 


MY, 7  )  =  -V 


»-  =  c.osrc 
v  c  =  i .  r  c 


vi=nF>p(-r  ,rcnr ) 
y?-rfxp(-o  .  lcno ) 


r      n    fl    Ti-P    tct/u    nlFREP    cf    stfps 

C        T     I  <=     TH^    CLRPFNTTINF      


AI.Pt-A  IS 

N=2 


7FF  STFP  rt-ANGING  FAC7CR 


t  =  ?.dc<-f 
ALP»-A  =  o.ino 


7 LPT  t-Nl  ,i-\2  ,  h  N3  PE  7FPFF  CCNSFCU7IVE  STEPS  CVER  THE  TIME 

f    INTFRVAL  TN1»TN2,TN3  ^_ 


LFT     FFl,FP2tFP* 
HM=O.C=CC 


FE     TFE    CCEFFICIENT    CF    A-P     FPECICTCR 


l-fOrCC^C 

HN^=0.05CO 


C       LET    YN  ,VM  ,>N?,VM     PF     THE     VALUE    CF    Y    CVEP    THE     TI*E 
C        INTERVAL    TN  » TM  f  TN2  » TIS  3 


Y  M  =  Y  2 
YN"  =  Y  1 


YIS3-YC 

fCwFlTF     TF-E 


Pc-FT  ICT  FT  VALLF  CF  Y 


mc 


RPl  =  l.n:  +  (HNlM?.Cr<FM+6.Cr,*FN2  +  3.n*FN3))/(6.rr-'-FN?*(FN2  +  HN' 

pp7^-(hM^(2.CC:FM^Q.CC^FN2^a.nC^FN3)W(fc.CC^FN2^N3) 

pp^  =  ((-M^(2.C0J;l-M43.C0^FN2))/(6.nC*FN3*(FN2*  F  N  3  )) 


T=T+FN1 


TN  =  T 

TM  =  TN-HM 


7N9  =  TN  1-HN? 

tn^=t N2-fisq 


YN  =  YM+HM"(RFl-F(YM,TM»+nF?*F(YN2.TN2l  +  PP3*F(YN3,TN?)) 
LFT  SC0tSCltSC2»SC*  3F  7F":  CCEFFICIENTS  CF  THE  SIFF  CCPPECTCR 


nFL  =  (TN-TM)^(TM-T(\:,)->(TN2-TN?)«-((TNZ-TN3)^{TK-TNl)^ 

C (^tn-2-7N?-7N3)^TM-TN'»)»(T^TI^3>'M3*TN-2-'>TN1-TM-M  + 

~C  (  TM-TN?  r*l  TN2-TN3)'»-TM  +  TK2)  ) 

S  C  C  =  (  1  /  C  F  L  IMTM  -TM)»CTN2-TIS3)MTN-TN31»((TN?-TN?)* 


c  <TN-7Nn)"(7IS-7IS?)-<TM-TN?)*(TN-TN3)*(T!S-TM  )♦ 

_C (  TM-TM  ?  1  -»  (TIS7-TN?  )  M  7  N  2-7N  1  M 

S  f  1  =  (  I  /  C  E  L  »  -  (  T  N  -  T  M  )  ■  (  T  N  -  T  \  3  1 i  ■*  »  2  M  7  N  2  -  T  IS  3  )  * 
C  (  (  T\  2-TN:>)   •(?'■•  TN-7N2-TN?)-(TN-TN?)»»?) 


SC?  =  (l/r.  cl)-'(T\-TM)*(TN-TN?)**2*(TM-7N?)* 
f  (-(7M-TM)*  (?*TN-7M-7N3)-H7N-7N3)*»?) 


sc.  a=i  .  nc-cc  i-sr  ? 

C r F  P  P C  7     7  F  c     VALUE     CF 


Y    TUICF 


YN  =  sci'-YM  +  $c?y'VN?  +  $c3*YN3  +  FM*scc<F(YN»7M 

YN  =  SCl'YM-»<;C2;YN/  +  cC3;YN:>»HM*SCC*F(YNt7M 

PECCRF  7FF  TN7EGPA7ICN  FlSTfRY  EVERY  1C  S7EFS  ALCNC  7HF  WAY 

ip( (N/ir.nr-N/in > .f-T.c.rirc)  gc  7c  15c 


YFXAC7=CFXP(-7) 
P  P  T  N  7  c?c,7,YN,YFXAr.7 
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c?C       FCP»'AT(ri2.'5,2C?C.lQ) 


1  5  f       IF(T,G7.«,rCI    nc    TT    ?cc 
C        1MTI AU7E    ThF    HATA    NECESSARY    FCR    THE    NEXT     STPP 


hm=hn? 

hN?=HM 


FM=Al  PhA*HM 
Y>^  =  YN2 


V  N  *»  =  V  M 
YM=YN 


N  =  N  +  1 
CHANGE    ThF    STEP' SI7F 


ALPM  =  1  .PC/ALPHA 
GC     TC     KH 


20C       CC^TINLF 

STfP 

FNT 


8U 
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